» J 
{ THE QUARTERLY JOURNAL '¢ 


} MECHANICS AND 


| APPLIED 
| MATHEMATICS 


VOLUME VII PART 4 
DECEMBER 1954 


ENGINEER 
AAEERING 


OXFORD 


AT THE CLARENDON PRESS 
1954 


Price 15s. net 


R BRITAIN BY CHARLES BATEY AT THE UNIVERSITY PRESS, OXFORD 





THE QUARTERLY JOURNAL OF MECHANICS 
AND APPLIED MATHEMATICS 


Editorial Board 


S. GOLDSTEIN R. V. SOUTHWELL 
G.I. TAYLOR G. TEMPLE 
together with 

C. AITKEN H. JEFFREYS 
CHAPMAN J. E. LENNARD-JONES 

R. COLLAR M. J. LIGHTHILL 

G. COWLING G. C. McVITTIE 

G. DARWIN N. F. MOTT 
W. J. DUNCAN W. G. PENNEY 

E. GREEN A. G. PUGSLEY 

A. HALL L. ROSENHEAD 

R. 


A 
s 

A 
T 
Cc 
A 
D 


A. 
HARTREE O. G. SUTTON 
L. HOWARTH ALEXANDER THOM 
WILLIS JACKSON A. H. WILSON 

J. R. WOMERSLEY 


Executive Editors 
V. C. A. FERRARO D. M. A. LEGGETT 


THE QUARTERLY JOURNAL OF MECHANICS AND APPLIED MATHEMATICS is 
published at 15s. net for a single number with an annual subscription (for 
four numbers) of 50s. post free. 


NOTICE TO CONTRIBUTORS 


1. Communication. Papers should be communicated to one or other of the Executive 
Editors, by name, at King’s College, Strand, London, W.C. 2. 


2. Presentation. Manuscripts should preferably be typewritten, and each paper should 
be preceded by a summary not exceeding 300 words in length. References to literature 
should be given in standard order, author, title of journal, volume number, date, page. 
These should be placed at the end of the paper and arranged according to the order of 
reference in the paper. 


3. Diagrams. The number of diagrams should be kept to the minimum consistent with 
clarity. The lines of the figures should be drawn in ink either on draughtsman’s 

or on good quality white paper. Each individual line in the figure should bear reducing 
to one-half of the size of the original, and great care should be exercised to see that the 
lines are regular in thickness, especially where they meet. Letiering of the figure should 
be in pencil and should be sufficient to define clearly the lines and curves in it. The 
writing of formulae or of explanations on the diagram itself should be avoided. All 
explanations of symbols, etc., should be given in underline. Contributors should 
indicate on their manuscripts where figures should be inserted. 


4. Tables. Tables should preferably be arranged so that they can be printed with the 
columns parallel to the longer edge of the page. 


5. Notation. All single letters used to denote vectors in the manuscript should be 
marked by underlining with a wavy line. Scalar and vector products should be denoted 
by a.b and a » 6 respectively. Real and imaginary parts of complex quantities should 
be denoted by re and im respectively. 


6. Offprints. Authors of papers will be entitled to 25 free offprints. This number is 
available for sharing between authors of joint papers. 


7. All correspondence other than that dealing with contributions should be addressed 
to the Publisher: 


GEOFFREY CUMBERLEGE 
OXFORD UNIVERSITY PRESS 
AMEN HOUSE, LONDON, E.C. 4 





ptt 








BENDING OF A CURVED BAR IN ITS OWN PLANET 
By D. RADENKOVIC (Beograd, Yugoslavia) 
[Received 11 December 1953] 


SUMMARY 


The problem of the plane bending of curved bars is dealt with under the assump- 
tions that: (1) the influence of the deflexion on the conditions of equilibrium has 
to be taken into account; (2) the deflexion is, however, so small that the squares 
and the products of the displacements and their derivatives can be neglected. 
Linearizing the Kirchhoff—Clebsch equations of equilibrium and introducing a suitable 
fictitious loading, the basic integral equations of the problem are formulated. These 
equations permit a qualitative analysis of the problem and at the same time yield 
aready numerical analysis of any particular case. The application and the practical 
importance of the proposed theory is shown in a numerical example. 


1, Introduction 

In this paper we shall deal with the problem of the bending of bars from 
the point of view of a linearized theory founded upon the assumptions that 
(a) the influence of the deflexion on the conditions of equilibrium has to be 
taken into account; (b) the deflexion is, however, such that the squares and 
the products of the displacements and their derivatives can be neglected. 
We also assume that the material is perfectly elastic and that the central 
line of the bar lies in a plane which is also the plane of bending of the bar. 
The theory of bending and buckling of a straight bar resulting from the 
above suppositions was developed a long time ago. Yet no attempt has 
been made to formulate an analogous general theory{ for a curved bar 
which would permit a qualitative analysis of the problem and which would 


at the same time yield a ready numerical analysis of any particular case. 


2. Survey of necessary known formulae 
The well-known basic differential equations of the problem are: 
(1) Relations between the natural coordinates of the central line%p, s) in 
the original and deformed states expressed in terms of strains (e, 8): 
1 1 _ dp (1) 
P Po ds, 
(2) Relations between strains (e, 8) and stress resultants (M,N, 7): 
N dp M 


EA ds, El 


t This paper represents a condensed version of part of a D Sc.(Eng.) dissertation approved 
by the Technical University of Belgrade. 

t Some particular solutions can be found in the German technical literature; cf. for 
instance (1) and (2). 


ds ds,(1 +-¢€), 


(2) 


(Quart. Journ. Mech. and Applied Math., Vol. VII, Pt. 4 (1954)] 
5092.28 ce 
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(3) Relations between strains and displacements (wu. v): 


dv u du v 3 
«= ——-, = —-+-—. ; 
dsy Pro ds, Po : 
(4) Conditions of equilibrium (Kirchhoff—Clebsch equations): 
e+ +0 = 6) 
ds p 
qN TT. T =O}. (4) 
dsp 
dM 
‘Mg 0 
ds 


The meaning of the notations and sign-conventions are shown in Fig. 1; 
in the figure, as in the equations, the index zero refers to the original state. 





It is important to note that in the general case the components of the 
specific load (U,V) along the normal and the tangent to the central line 
may change during the deflexion. We shall consider here two cases of 
practical interest in which the intensity of the load remains unaltered: 
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(a) if the direction of the load turns together with the element of the bar 
(e.g. hydrostatic pressure), then U = U, and V = Xj; (6) if the loading 
keeps the original direction whilst the element turns (e.g. gravitational 
forces) then, as is readily seen, U = U,—V,B and V = K+ UP. 

Our problem is completely defined by the nine equations (1)-(4) with 
nine unknowns (s, p, «, 8, u, v, M, N, 7). Equations (1)-(3), being derived 
under the assumptions that the deflexions are small and that Hooke’s law 
is valid, are linear; whilst equations (4), expressing conditions of equili- 
brium of an element of the bar in the deformed state, are non-linear. 

Usually, in strength of materials and in structural analysis, the problem 
is linearized by writing the conditions of equilibrium with reference to the 
original position of the central line. We may express this by rewriting the 
equations (4) in the form 

dT, | No, U, : Q | 
ds Pg 
dN, Ty 


Is T Ny 0}. (45) 
ASo Po 
dM, 
a 0 
ds, 0 


Equations (4,) together with (2) and (3), equations (1) being now super- 
fluous, form a complete set which is the basis of conventional structural 
analysis. In the following we shall refer to a solution of this set of equations 
as ‘the solution for a rigid bar’, although in the analysis of statically 
indeterminate systems the displacements appearing in the end conditions 
must be taken into account in order to find the stress-resultants. The 
strains and stress-resultants in the rigid bar will be denoted by the index 
zero, aS has been done in the equations (4,). 

Actually, in structural analysis the set of equations (2), (3), (4)) is never 
used, because the integrals of these equations can be immediately written 
down in any particular case from purely statical considerations. Thus the 
conditions of equilibrium can always be expressed in terms of stress- 
resultants (involving no derivatives), and the displacement, say &(s), of a 
point S in a certain direction is, by virtue of the principle of least work, 
always given by 


l l 
* M(r)M,(r, 8) r N(r)N,(r, 8) 
E(S —~ dr + - fia Ir, ‘ 
()= | re + | ae ™ (5) 
0 0 


where M(r) is the actual bending moment at any point r; M,(r,s) is the 
bending moment at 7 produced by a unit load at s acting in the direction 
of the required displacement; NV and N, have similar meanings. If € denotes 
the rotation of the cross-section, then M, is produced by a unit moment at s. 
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Moreover, as the problem is linear with respect to the loads, any unknown 
quantity (e.g. M, B, », €) denoted by k(s), produced by a given load q(t), 
can be found by means of an influence line: 

l 
k(s) = | k(s, t)q(t) dt, (6) 

0 
where k(s,t) is the value of k at the point s produced by a unit load at t 
acting in the direction of the given load q(t) at that point. Now, k(s,t) is 
the Green’s resolvent of the corresponding differential problem defined 
by (2), (3), (49). But, as is seen from the above, this resolvent can always 
be found from statical considerations without reference to the underlying 

differential problem. 


3. Linearization of the conditions of equilibrium 

After this necessary survey of the well-known facts we may pass to our 
actual problem. 

We now wish to express the influence of the defiexion on the conditions 
of equilibrium by equations (4), but, because the deflexion is assumed 
to be small, we can linearize these equations by neglecting the terms of 
second or higher order. Thus we can neglect the influence of the deflexion 
due to the extension of the central line on the conditions of equilibrium, 
i.e. we can put ds = ds,, which of course does not imply neglecting the 
effect of the extension on the end conditions in a statically indeterminate 
structure. Also, by virtue of the second equation in (1), we have 

N y( * )- 

p Py 189 
But as the actual value of N differs very little from the normal force JN, in 
a rigid bar, say N = N,+N, and as the product N(dB/ds,) can be neglected 
as a small quantity of higher order, then we can put 

N WN, dp 

+-% 

Pp Po d8q 

"om 
and similarly : = : 4 7B 

, p Po dso 

With this we can rewrite the equations (4) in the form 


aT Ny Bp _o) 


J = 0 
dsy Po ° ds, 
aN T 1B : , 
ie |. ee oe (4') 
dsy Po ds, 
dM _T=0 


d So 
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These equations are linear with respect to the unknowns MV, N, T’, B, but 

they are not linear with respect to the external loading. We have now to 
take into account the behaviour of the external loads. 

(a) If the direction of the load turns with the element of the bar, then 

U = U, and V = J, and as d8/ds, = —M/EI the equations (4’) become 


0 0? 
dT N 1 (No ay\ 0 


ds,' pp. ° \#I J 
wT .,{@ea 
ae 0M\=0>}. 4a 
an a ee te 
dM _T=0 
dso 


(b) If the loading keeps the original direction, then U = U,—\B and 
y=V,+U,8. Multiplying the equations (49) by 8, we find the values of 


J,8 and U,8, which, after some rearrangements, enables us to write the 


equations (4°) in the form 


aT + u,4{£ p)——28! = 01 
dsy Po \dsq Po | 
aN | | V | d (T', B) n No a Qs. (4 b) 
ds, Po \dsy Po 
dM _. 0 
ds, 


From the linearized equations of equilibrium (4a) and (4b) we could 
easily deduce the already known differential equations of bending and 
buckling of a straight beam. The differential equations for special cases 
of buckling of arches investigated by Lockshin (3) and Chwalla (4) can 
also be obtained from them in a simple and general manner. However, 
this would lead us to differential equations of the sixth order with variable 
coefficients, which, especially in view of the number of end conditions, 


would be too complicated to be dealt with. 


4. Basic integral equations of the problem 

Now it can be easily seen that the equations (4a) and (4b) are identical 
with equations (4,) if the bracketed terms are considered as a fictitious 
load (U,, V, and U,, V, respectively).+ Thus it is quite clear that all the 
integrals occurring in the rigid-bar analysis remain valid for a bar under 
the actual and fictitious loads, the set of differential equations of the 
problem being the same in both cases. In particular, we can use Green’s 


+ It may be noted that Timoshenko (5) has already used the notion of the fictitious 


ud, but in a very special case and with quite different aim. 
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resolvent of the problem, as indicated in equation (6), which will bring us 
to the non-homogeneous Fredholm’s integral equations of the second kind 
as the basic equations of the problem. These equations will permit us to 
formulate the problem of a plane bending of curved bars in a way which 


—N, M 
El 





is most general and at the same time best fitted for both qualitative and 
numerical analysis. 
Let us consider the two cases of the loading mentioned previously: 


(a) If the direction of the load turns, we shall write the basic equation 
in terms of bending moments. From the equations (4a) and Fig. 2 it is 
seen that in this case the fictitious load q,(¢) is at all points perpendicular 
to the stress-resultant F, for the rigid bar and has the intensity F,(M/E1). 

Then, writing the equation (6) for the moments and noting that 


I 
' M(s, t)qo(t) dt M,(s), we get 


0 





l 


M(s) = Ma(s)+a { M(s,t) aa M(t) dt, (7) 


0 

where A is a parameter varying with the intensity of the loading, and 
M(s,t) is the bending moment at s produced by a unit load at ¢ acting in 
the direction of the fictitious load q,(¢). It will be noted that the kernel of 
the integral equation (7) is a product of two factors: ‘influence line’ M(s, t) 
—main part of the kernel non-symmetrical in the general case—and 
F,(t)/ EI (t); both of them, being defined for a rigid bar, can always be found 
by usual methods of structural analysis. 

It is worth noting that in numerical computations it may be of help if 
the fictitious load is resolved in two fixed directions. say x, y, so that instead 
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A CURVED BAR IN ITS OWN PLANE 





BENDING OF 





of M(s,t) depending at every point on the direction of F(t), the influence 
lines ,(s,t) and M,(s,t), referring at all points to these fixed directions, 
are used. 

We may mention that Smirnov (6) formulated in terms of deflexions 
similar integral equations for the problem of buckling of arches. 


(b) If the loading keeps the original direction, we are bound to write the 
basic equation in terms of the rotation of the cross-section 8. But in this 
ease the direction of the fictitious load is not defined by the direction of 
the stress-resultant F, in such a simple manner as before, nor, as we shall 
see, is it important here. To begin with, we shall deal with each com- 
ponent (U;, V,) separately. Then, denoting by £,,(s,¢) and £,(s,t) rotations 
at s due to a unit force at ¢ acting in the directions of the normal and the 
tangent, respectively, of the centre line at t, and keeping in mind that 

; 
B.(s) = | B(s, t)q,(t) dt, we get, in virtue of (6): 


0 
i 


B(s) = B,(s) ra | 


l 
Ba(8s1) 5(NoB)—Bels.t) 5, (Top) dt — 


Ty(t) , N,(t) 
—A , + B,,(8, £) —— d 
fe Ex 4 Pol) Bela, ae . 


0 
By(s)+A(1,+1,). (8) 


In order to transform this integro-differential equation into an integral 





equation, we integrate J, by parts: 
l 


I, = [B.(s, ()No(t)B(t)—8,(s, t)To(t)B(t)], 4+ 
l 
, t) 
[ | — Pal nye) + Pe) aye] (0) dt. (9) 
ot ot 
0 
Now, at supported ends f(s, 0) = B(s,l) = 0, and for free ends N, = 7, = 0. 
(If there is a force acting at a free end, it can always be treated as an external 
force acting very near the supported end.) 
Thus, on (9) in (8), we have 


a) Cc B, AS; t y ( v &,t 
Bls) = Bo sa f ee '(—Ny)4 = 7) 


0 
éy N, 
a B,(8, t) z o_ B.(8, t) re : 
Po Po 
The equation (10) can be simplified as follows. By Maxwell’s reciprocal 
theorem, the rotation at s due to a unit force normal to the central line at t 


|e dt. (10) 
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must be equal to the displacement normal to the central line at ¢ due toa 
unit moment in s; thus B,(8,t) = wat, 8), 


OB, (s,t 0U,,(t, 8 
so that Bus, t) — Cult, 2. (11) 
ot ot 
Now, by virtue of the second of the equations (3), 
ou (t, 8) ~.. 
mart = By (t, 8) Vy (t, 8). (12) 


0 


Inserting (12) in (11) and applying again Maxwell’s reciprocal theorem, we 


get Op,,(s, t) 
mane = By(s, t) )— Bel 8, t). (13) 
In the same manner we obtain 
OB..(s, 0v,,(E, « ] i 
Br(s, t) _— evult, 8) — €y,(t, 8) + is -Ux,(t, 8) — ——— B,(8, t) (14) 
ot ot Po Polt) 


aS €y,(¢, 8) can always be neglected, being a small quantity of higher order. 
By the use of (13) and (14) the equation (10) reduces to 


B(s) = Bo(s)+A { Bu(s.t)[—Ny(O]B() at. (15) 


Here again we have an integral equation of the same kind as (7). The 
kernel of the equation (15) is also given as a product of two factors, but it 
can be shown that here the main part of the kernel f,,(s, t) is always sym- 
metrical and positive definite. Since f,,(s,t) is the rotation of the cross- 
section s due to a unit moment at t, by virtue of the equation (5) it always 
yields the integral representation 


By,(s, t) ate Nd Anes a EY (16) 


I 
_ fF Mr, 8) M(r, t) 
‘ Eli(r) 


. 


0 
where IM(7,7) is the bending moment at r due to a unit moment at 7. The 
above properties of 8,,(s,t) are explicitly seen from (16). It is interesting 
to note that, for the case of a straight bar, Trefftz (7) developed an equation 
identical with (15) by a different procedure which would hardly permit the 
generalization obtained above. 

It may be mentioned here that all the above analysis is valid only for 
distributed loading. If the bar is loaded by concentrated forces, or if the 
end conditions are given with respect to the stress-resultant Fj, fictitious 
concentrated loads may appear in the computations. This would change 
the form of the kernels in the general case, but of course the principles of 
the analysis would remain the same. Because of lack of space, we shall 
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omit here this interesting but rather lengthy analysis, which will be pub- 
lished later with the full text of the author’s dissertation. 


5, General analysis and numerical solution of the problem 

From the theory of integral equations, it is well known that an equation 
of the type (7) or (15) possesses a unique solution, except for special values, 
so-called eigenvalues, of the parameter A. In order to find these values the 
corresponding homogeneous equation, obtained by putting M,(s) and £,(s) 
equal to zero in (7) and (15) respectively, has to be solved. The homo- 


geneous equation corresponding to (7) is 


F(t) 


BQ uo 4 (17) 


M(s) =A | M(s,t) 
0 
and the homogeneous equation corresponding to (15) after some trans- 


formations can be put in the form 


l 
M(s) =d | G(s,t)M(t) dt, (18) 
0 
l 


[ Ms, r) M(t, 1) 


enn Welt) Or. 19 
El(r) o(”) dr (19) 


where G(s, t) 

0 
Now the well-known method of Vianello for the computation of critical 
loads is actually generalized by the equations (17) and (18) to the case of 
curved bars. The solution of the non-homogeneous equation given by 
Neumann’s series, which converges for |A| < |A,!, A, being the least eigen- 
value, when applied to our problem, represents another generalization of 
a method commonly used in the computation of straight beams. 

In connexion with the problem of convergence of Neumann’s series, but 
interesting from many other points of view, is the question of the existence 
of the critical load, i.e. of the reality of eigenvalues. 

Now, as we have shown that the main part of the kernel of the equation 
15) is always symmetrical and positive definite, it follows from a well- 
known theorem in the theory of integral equations that, if the loading 
keeps its original direction, then real and only real eigenvalues always exist. 
With this, Trefftz’s theorem on the existence of critical loading for straight 
beams (7) is generalized and at the same time made more precise, its validity 
being restricted to this case of loading. 

By virtue of a theorem due to J. Marty (8) it is possible in some cases 
to symmetrize the kernel of the equation (7) as well, but it can be shown 
that, if the load turns with the element of the bar, it may happen that the problem 
has no eigenvalues, i.e. that no critical load whatever exists. This happens 
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always if the loading acts only on that part of a bar between a free end and 
an immovable support which is next to the free end (e.g. cantilever beam), 

The analysis of the variation of eigenvalues with the small change of the 
kernel would lead to other qualitative conclusions, which we have not 
considered so far. 


qi: 

























































































Fic. 3. 1 = 212-:00m.; f = 21-25m.; g, = 8-80t./m.; g, = 13-00 t./m.; 
I, = Icosdé = 0:-460m.4; F = 0-319 m.2; E = 2-1 10% t./m.? 


There remains the question of the actual numerical analysis of any 
particular problem. It can be dealt with in the following well-known way. 
As we can approximately evaluate an integral using a finite sum, any 


integral equation of the form 
l 


y(x) = f(x)+A | K(x, z)y(z) dz (20) 
0 
can be replaced by a system of algebraic linear equations 
y = f+A, KAy, (21) 
where K is a matrix corresponding to the kernel of the equation (20) and 
A is a matrix depending on the adopted method of numerical integration. 
If any classical method of numerical integration is adopted (e.g. trapezoidal, 
Simpson’s, Gauss’s, Tschebysheff’s) the matrix A will be diagonal. In our 
case the method of ‘elastic weights’ known from structural analysis, in- 
volving a continuant matrix A, yields the best results. Further details and 
numerical examples will be published in the author’s dissertation. 


6. A numerical example 

We shall illustrate the application of the proposed method by the 
example of the analysis of a bridge arch shown in Fig. 3.+ 

The load is supposed to turn during the deformation, as in this case the 


+ In order to compare the results the same numerical data are taken as in Fritz’s book 
((1), p. 34). Fritz treated the problem by a method usual in German technical literature 
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basic equation (7) of the problem is given in terms of the required bending 
moments.t We shall write this equation in a slightly different form because, 
as has been mentioned already, it is easier to evaluate the influence-lines 
with respect to two fixed directions. So we have 


l 
7 F,,(t) F,,(t) 
M(s M,(s)-+-A M,,(s, t)—¢— — M,(s, t)—°=— | M(t) dt, 
(8) = Mla)-+0 | | Ae, et) — Mylo, 0) 52 |e (a) 
or if we take for A the constant factor —F,,/Hj, i.e. putA = H/EI,, 
I 
M(s) = M,(s)+A [ x s,t)—* M(t) dt, b 
(s) \(8) ( Tt (t) ¢ (6) 
0 
where K(s, t) M,(s, t)+-M,(s, t) a (c) 


Here 1/,(s,t) is the moment at s due to a unit load at ¢ acting in the 
direction of the x-axis, M/,(s,t) has a similar meaning, and /), and Fj, are 
the components of the stress-resultant K, in the directions of 2 and y 
respectively. For the numerical analysis we replace the equation (a) by 
the system of linear algebraic equations 


M = M,-+A, KAM. (d) 


For this purpose the arch is divided into eight parts (with equal horizontal 
projections), so that seven unknown bending moments appear in the 
equations (d) and the matrix K is of the seventh order. The free terms M,, 
and the elements of the matrix K,,,, = K(s,,,8,) are found by the usual 
methods of structural analysis which supply us with the following data: 


Stress-resultants for the given loading 


























| 
I I 7 5 4 5 6 7 
Uy 2371 222 2554 365 1870 — 2677 2053 
| 
Fy g22°2 577°7 233°2 |—111°-3. |—344°5) |—577°7 —810°9 
l — 2864°5 
Influence lines for the vertical loading M,(8,,8,) 
—- ee | 
I 2 3 4 | 5 6 7 
16-198 | 7042 |—o-108 |—4:757 | 6-733 |—6:208 |—3-677 
7°893 | 17°75 4°547 4°370 |—8-703 |—8-749 |—5°357 
rcSc 5°627 13°965 1°163 |—5-910 |—7°623 | — 5040 
726 2°831 1°647 11°84! 1°647 |—2°831 | 2-726 





T Actually here the load (due to gravitation) keeps the original direction, but in fact it 
an be show n 


that the behaviour of the loading seldom influences the results appreciably. 
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Influence lines for the horizontal loading M,(s,,,, 8,,) 











n | | | 
= 
m | I 2 | 3 | 4 | 5 6 7 
\— --——— —— —- -— - - Oe 
I 5°662 3°346 | 2°247 1°992 | 2°069 1°966 | 1311 
- 2°734 5°167 | 3°142 2°656 2°836 | 2-802 | tors 
3 o'sir 1°478 | 2°682 I'992 | 2°299 | 2-506 1°813 
4 —1°004 |—1'079 |—0-°460 0000 | 0-460 | 1:°079 1°004 





With these data from the usual structural analysis we can pass to the 
actual analysis of bending taking into account the deflexion.+ 
First we evaluate the elements 


Ran a M, (8m; 8), ) + M,(8; 8p ) vou) : 


which can be immediately entered in the following table by means of a 
computing machine. 


The matrix K 














n | 
a 

m is I oa 3 4 5 6 7 
I | 18021 | 7°717 0075 |—4°834 |—6-982 |—6-604 |—4-048 

| oe | - » oo — | 2 

2 | 8773 | 18793 4°803 4473 |—9°044 |—9°314 |—5*899 
3 |} 1*750 5°925 14183 | 1°086 — 6°186 | - 8128 | —5°553 
4 | — 3°049 -3°049 1°610 11°841 1°592 |—3'049 |—3:010 
5 |—5°624 |—8-128 |—6-097 1°240 | 14°288 5°925 1°732 
6 |—5°974 |—9°314 |—8°934 [—4°207 | 4°915 | 18°793 | 8-667 
7 | —4°099 |~6-604 — 6-901 —4°680 | o162 | 7°717 17‘801 








For the integration we use the ‘method of elastic weights’ and assume 
that the influence line between any three consecutive points is a parabola. 
Then, because in this particular example J(t) = J,sec¢ and at the same 
time dt = dxsec ¢, we have (cf. (6)) in the equation (c) 


5Al) 5ALH 


A = — - — 
a 6 6EI, 
10 O-]1 
0-1 41:0 O-1 
and A 0-1 1:0 O71 
0-1 1:0 


We rewrite the equations (d) in the form 
a 
Xo 


+ It may be stressed here that the following computations are given in only a very 


(- a kA)M M,. (e) 
Ao 


slightly abridged form, so that the total amount of work required to obtain the corrected 
values of bending moments can be readily appreciated from the following tables. 
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Taking into account that in this instance 1/A, = 152-708, we can imme- 
diately write down this system in the explicit form as follows: 











ve M. VM, M, M, M, = M,/r, 

I 527 4 5°525 8125 7°77 4°708 362222: 
0 32°555 6°235 4°909 10°422 10°808 6°830 492023°9 
7°518 137°824 1°885 | 6:8g0 9°302 6°366 390015*2 

"193 2°48 140°547 2°471 3°191 3°315 55738°3 

437 *300 786 2°059 | 137°704 7°527 2°325 | —255503°2 
10°805 10°292 4°669 6°368 | 132°557 | —10°546 | —408798-°2 

7 7°794 5029 5°354 0°466 9°513 | 134°135 | —313508°7 





As the diagonal terms are far greater than all the others this system is 
easily solved by means of successive approximations. 
The results are given in the first row of the following table 








2 3 4 5 6 7 
| e eq. (e) 2604 4932 3707 399 3051 — 439° — 331 
I 5 4956 3881 281 — 2837 4173 — 3150 
[ l 3222 2554 365 1870 2677 — 2053 

I o°5 3°0 27°7 7° 52 4°7 

34°7 32°2 6:2 38°7 39°0 37°9 





The bending moments computed by Fritz and the bending moments 
obtained by the usual structural analysis} are given in the second and the 
third row of the table respectively. 

The fairly good agreement between our results and those of Fritz is due 
to the fact that we have chosen an example in which all the assumptions of 
Fritz’s method about the form of the centre-line and the variation of the 
cross-section are fulfilled. For any other form of arch the method usual 
in German technical literature would give less reliable results, whilst the 
method described above leads to no reduction in accuracy and no increase 
in the work 

[t will be seen that the actual bending moments are 30-40 per cent. greater 
than those given by the usual structural analysis. Because the given data 
could apply to an actual steel bridge with wide span, we can easily see the 
practical importance of a reliable theory of bending of arches which takes 
into account the influence of the deflexion on the conditions of equilibrium. 
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t Which, as is already mentioned, implies neglecting the influences of deflexions and 


aking (45) as the equations of « yuilibrium. 





BENDING OF A CURVED BAR IN ITS OWN PLANE 


REFERENCES 

1. B. Frrrz, Theorie und Berechnung vollwandiger Bogentrdger (Berlin, 1934), 

2. F. Discuincer, ‘Elastische und plastiseche Verformungen der Eisenbetontrag. 
werke und insbesondere der Bogenbriicken’, Bawingenieur (1939). 

3. A. S. Locxsuin, ‘On the buckling of a rod naturally curved’, Philosophical 
Magazine (1932). Cf. also: 

3a. A. N. Dinnik, Ustoichivost’ arok (Stability of Arches) (Moscau, 1947). 

4. E. Cuwatta, and C. F. KortiBrunner, Uber das Ausknicken symmetrischer 
Bogentrager unter symmetrisch verteilten Belastungen (Stahlban, 1937). 

5. S. P. TrmosHEnkKo, Theory of Elastic Stability (New York, 1936). 

6. A. F. Smirnov, Staticheskaia i dinamicheskaia ustoichivost’ sooruzhenii (The Statical 
and Dynamical Stability of Structures) (Moscau, 1947). 

7. E. TreFrtz, ‘Allgemeine Theorie der Knickung der geraden Stabes’, Z. f. angew. 
Math. u. Mechanik (1923). 

8. J. Marty, Three notes in Comptes rendus, 150 (1910). 











DY 


Th 
Maxv 
with « 
The ° 
modu 
relax 

Ex 
obtai 

Th 
stres: 

(a. 
than 
belo 

(b 
solid 
struc 
dyns 

(c 
ina 

T 
shoy 
acr 

dyn 


" 
TH 
boc 





mMtrag. 


ophical 


trischer 


Stat 1 cal 


angew, 











DYNAMIC BEHAVIOUR OF LINEAR RHEOLOGICAL 
BODIES UNDER PERIODIC STRESSES 





By 8. IRMAY?* (Israel Institute of Technology, Haifa) 
Received 15 October 1953] 


SUMMARY 

The simpler linear rheological bodies—such as the elastic and Kelvin’s solids, 
Maxwell’s and Newton’s viscous fluids—obey a generalized linear differential equation 
with constant coefficients between the strain and stress deviators and their time rates. 
The various coefficients are correlated to the asymptotic rigidity or static elastic 
modulus G, elastic firmness or dynamic modulus H, solid viscosity pu, times of 
relaxation R and of lagging or retardation L, and the endosity 7. 

Expressing the time and stress in non-dimensional forms, a universal equation is 
btained, dependent on a single non-dimensional parameter, the ‘time factor’ 
+= R/L = G/H. This defines a principle of similitude for all bodies of identical 7. 

The mechanical and thermodynamic study of periodic, impulsive, and transient 
stresses leads towards a new classification of the linear bodies, based on the value of r: 

a) The exothermal or dissipative bodies (r < 1) are less rigid than firm, less strained 
than elastic solids, dissipate energy in a cycle. They warm up adiabatically. Here 
elong also Maxwell’s elastico-viscous fluids and Kelvin’s firmo-viscous solids (7 = 0). 

b) The endous solids (+ |) are more rigid than firm, more strained than elastic 
solids, may show a ‘negative’ dissipation in a cycle, which is compensated by internal 
structural changes. This is not in contradiction with the second law of thermo- 
lynamics and may explain the fatigue of such solids. 

c) The homothermal solids (r 1) are as rigid as firm, show no apparent dissipation 
nacycle. A particular case is that of elastic solids. 

The stress-strain diagram under periodic stresses tends towards an ellipse, which 
shows accommodation of the body. The smaller axis passes through a maximum for 
critical frequency, and the body has two independent elastic moduli—a static and a 
lynamic one. These facts were experimentally confirmed for certain plastics. 


1. Introduction 

Tue rheological equation of a generalized homogeneous isotropic linear 
body with constant scalar rheological coefficients§ a,, dy, dg, a4, relating 
the four deviators of stress p, strain e, and their time rates » = ép/ét and 


€= Ge/ct (flow deviator), is 


a, p+a.p = dge+ay,é. (1) 
A brief summary of this paper was presented orally by Professor M. Reiner at the 
August 1952 Congress of Theoretical and Applied Mechanics, Istanbul. 
{ These terms are defined in the Introduction and justified later. 
§ Hohenemser and Prager ((1), compare Reiner (2)) add here a constant deviator c, 
hich means that, without stress and whatever the past history of the body, there is always 
same f strain « Replacing « by «—e,, then writing « for e—e , we get (1). The 
stant term of the plastic Bingham body p c+ pe is not a truly linear term, as the 
uation is valid only for p 


(Quart. Journ. Mech. and Applied Math., Vol. VII, Pt. 4 (1954)] 
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Dividing by any one of the non-vanishing coefficients, we obtain four 
parallel equations, each characterized by three independent coefficients 


p+Rp = Get+pe = G(e+ Le) = p(e/L+e) 
p/R+p = e/nt+He = H(e/L+€) 
p/G+np = «e+Le 
p/pt+p/H = e/L+e J 
Six ratios and a non-dimensional parameter enter here, the names of which 
will be justified later on, viz. 


G =a,/a, (dimensions F L-*) is the static modulus of elasticity | 
or asymptotic rigidity 

H =a,/a, (dimensions F L-*) is the dynamic modulus of elastic 
firmness 


R= a,/a, (dimension 7’) is the time of relaxation 
L = a,/a, (dimension 7’) is the time of lagging or retardation 
fe =a,/a, (dimensions FL-*7) is the solid viscosity 


7 = 4a,/a, (dimensions F-'L*7') may be called endosity 


T = G,4,/a,a, (no dimensions) may be called the time factor | 
As only three coefficients are independent, we have the relations 

G = R/n = p/L; R= p/H = Gn; p= GL= RBH | 

H = L/q = p/B; L = p/G = Hn; n= R/G=L/H +. (4) 

7= R/L=G/H J 


2. Principle of similitude 
Equation (1) is symmetrical with respect to the origin of the stress-strain 
diagram, linear, i.e. independent of the scale chosen, and additive. 
When L is taken as unit of time and G as unit of stress, we can express 
the time, stress, etc., in non-dimensional forms, which may be called the 
numerical time 6, stress P, flow «’, and stress time rate P’, by 


@6=t/L; P=p/G; &=0/00=—Lé; PP’ = @P/200= pL/G. (5) 
Equation (1) becomes a universal or numerical equation with a single 
parameter 7 (time factor) 

P+rP’ e+e’. (6) 
The value of 7 is characteristic of a group of linear bodies behaving similarly. 
This defines a principle of similitude. 
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3, Stress-strain diagrams 
If, starting from rest (« = 0; p = 0 or P = 0), we apply a stress p(t) or 
the numerical P(@), the solution of (5) is 


« = p/H+(1/G—1/H)exp(—t/L) [ p(t’ exp(t’/L) dt’ 
t’=0 
6 
tP+(1—s)e® | P(6')e® dé’. (7) 


6’ =0 
































IW 
} w=0 
iy, 
\¥ r=() 
| 
1 
WwW | @W 
w T=2>1 
\¥ 
T= 
1 - T 
WwW WwW 
Fic. 1. Amplitude factor WJ. 
An important case is that of periodical stress 
p=p,sinft or P= Psinw, (8) 
where f (= 27/7’) is the frequency, w (= Lf = 27L/T) is the numerical 


lrequency, 7' is the period, p, is the stress amplitude, P, = p,/G@ is a 


humerical stress amplitude. 


pd 
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The strain (7) becomes 
€ = e,[sin(w0+ «)—e~® sin a] 
= P,(1+w?)-1[(1-+7w?)sin w§—w(1—7)(cos w—e-)], (9) 
a= Pb b= V{(lt+re%(1 +04}; 


tana = —(1—r)a(1 +72) 


where 

(10) 
The exponential term vanishes for ¢ (or 0) > 00 and we get finally a pure 
sinusoidal strain with amplitude «, = P,y% and phase difference «. 


€ = €,sin(w+ a) = e,sin(ft+-a). (11) 


Stress Y~p 
+ 1+ 





‘8+ T=3 
‘65 w=!) 
a tana=0°5 
3 + 
pY™ . Strain 





1-864 -2 O42 47-6 8410 &, 





Fic. 2. Stress-strain diagram for periodic stress. 


Both a and y depend on the body (r) and numerical frequency (w) according 
to Fig. 2. ¢(7,w) is represented on Fig. 1. 


For 7 = 0, x = (1+w?)-4. 
For 721, pel. 
For rt > 00, yb > tw(1+w?)-4. 


At low frequencies (w — 0), 4 — 1; at high frequencies (w — 0), b>. 

Fig. 2 represents the numerical stress-strain diagram. We see that the 
curve tends asymptotically towards an ellipse; there is accommodation of 
matter. The equation of the ellipse is obtained by eliminating the time! 
(or @) between (8) and (11), 


(P/P,)?—2(P/P,)(€/e,)cos a+ (e €,)? = sin’. (12) 


The axes of the ellipse bisect the strain and stress axes. The units of 
numerical strain and stress are «, and P, = p,/G respectively. 
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Ifa is the half-axis, positive bisector, and b is the other half-axis of the 


ellipse, we have: 


For t= 1, b = 0, a straight line 


For + l,a V2 cos $a; b = vV2sinja <a>. (13) 
For 7 < 1, a= v2sin ja; b = v2cosfa >a 


The smaller half-axis varies between 0 (for w= 0 and wo) and 
vr—1|(7-+1)-* (for w 7~'). This fact has been observed experimentally 
by L’Hermite (3) for certain plastics. 

In the dimensional stress-strain (p,«) diagram we have also an ellipse, 
the a-axis of which is inclined at an angle y to the e-axis, with 


tany = G/b = G(1+w?)*(1+7?w?)-! = & following 7 = 1. (14) 


+ 


For w > 0, tany > G ) 
(15) 
For w > 0, tany>G/r = H 


We see that G and H are both elastic moduli: G is a static modulus, valid 
for constant or low-frequency stresses; H is a dynamic modulus, valid for 
high-frequency stresses. 

This may serve as basis for the experimental determination of G and H 
or 7) of a linear body. The (p, e) diagrams for a very slow and a very rapid 
periodic stress give two ellipses, the slopes of whose axes give G and H 
by (15). 

When the periodic stress does not change sign and varies between 0 
and P,, we obtain again for f+ oo an asymptotic ellipse similar to that 
of Fig. 2, only P has to be replaced by (P—0-5P,) and ¢« by (e—0-5P,)/p. 


4. Transient stresses 

Transient stresses occur whenever stresses are applied or removed 
suddenly. 

(a2) When the time of application is very short and the stress intensity 
is very large, we may approximate the stress by Dirac’s impulsive function 
0(#) defined by 


0 for 0~0 
0(0) 4 (16) 
0 for é=0 


| 3(@) dé = 1. (17) 
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It is represented on Fig. 3a, 


P = 6(0); € = ete, = 78(8)+(1—r)e~. (18) 


(4) 
o (a) 





a? (b) 








Fic. 3. Dirac’s 6(@) and 
Heaviside’s y(@) functions. 


As for elastic (Hookean) solids, ¢,, = 5(@), the elastic strain ¢, $ €,, following 
721. For >, «,>0. 
(b) When a constant unit stress is suddenly applied and maintained 
(Fig. 3b), it may be represented by Heaviside’s unit function Y (6); 
8 
P=Y(#)= | 8(x) da; € = e+e, = Y(0)—Y(6)(1—r)e~*, (19) 
whereas for elastic solids «,, = Y(@). Initially « = rY(@) =, following 
7 = 1; finally « > 1 for 6+. 

(c) When a constant stress is suddenly applied (@ = 0) and removed 

(6 = 0,) (Fig. 4), 
P = Y(6)—Y(0—86,); «(0<0)=0; 
«(0 <0< 8,) = 1—(1—r)e-9; e (9 > 0,) = (1—r)(e? —1)e~*. 

(d) When periodical impulsive stresses 5(n0,) are applied at intervals 
6, = t,/L we get for @ > 00 a finite strain «, = (1—e~®)/(e%—1) + (e%—1)". 
At high frequency (@, small) «, > 1/0, = p, L/Gt,. If the stresses are 
interrupted at 6 = 6,; « = e~%(e%—1)/(e%—1) + 0 for 0 > oo and there is 
no permanent strain. 


(20) 


(e) When a constant unit stress is suddenly and periodically applied 
(9 = n@,) and removed (6 = n6,+6,), we get after a long time: 
Shortly before 6 = n6,, € = (1—r)(e*—1)/(e%—1), 
and shortly afterwards e = 1—(1—7r)(e%—¢%)/(e%_ 1). 


Shortly before 6 = @,+-6,, « 1—(1—r)(e%-1_ ])/(e%_ 1), 


and shortly afterwards € = (l—7)(1—e—")/(1—e—%)., 
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Rach time there is a step or jump Ae = +. This suggests an experimental 
method for determining the time factor 7 of a body: a constant stress is 
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Fic. 4. Constant stress suddenly applied and removed. 


periodically and suddenly applied and removed. The jump in the strain-— 
time diagram represents r. The nearer 7 approaches 1, the more the body 


resembles an elastic solid. 


5. The rheological coefficients 

In section 3 it has been shown by the body’s behaviour under periodic 
stresses that G represents a static and H a dynamic modulus of elasticity. 
These two moduli are independent of each other and their ratio G/H = + 
defines the characteristic ‘time factor’ of the body. 

Another interpretation of these moduli may be given. A linear body with 
a, = 0 defines Maxwell’s viscous fluid, a generalization of the ordinary 
Newtonian viscous fluid. By (3) G = 0,7 = 0. When a constant stress p, 
issuddenly applied (at ¢ = 0), 


‘ p,/H+-p,t/p. 


When it is suddenly removed (at t = ¢,), 


© py ty/p. 
Here p,/H disappears with the stress, so it is the elastic strain. Therefore 
H may also be called the firmness modulus; p, t,/u is the permanent strain, 
ut, replacing G (or H). Therefore, though without elastic rigidity (@ = 0), 
the Maxwell fluid has elastic firmness: it is a firmo-viscous fluid. The term 
elastic’ should preferably be conserved for elastic rigidity alone. This 
distinction is usually not recognized in rheological literature, though 
L’Hermite (3) mentions it. 

















406 S. IRMAY 


+0 
On the other hand, an impulsive stress of impulse value J = | pd 
t=—o 


gives a permanent strain « = J/u. Hence » = I/e, which corresponds to 
G = p/e of the elastic solids. The viscosity p is related to stress impulses 
as the rigidity modulus G is related to stresses. This provides for a new 
definition of the viscosity, valid for all Maxwell’s fluids. 


6. Thermodynamic considerations 
The power w expended by the deviator stresses in the deformation of 
the unit volume of the body is 


w= pe. (21) 
The power expended is stored by the elastic part of the body, and dissipated, 
i.e. transformed into heat, by its viscous part. In an isothermal process 
this heat must be extracted from the body during an exothermal process, 
or introduced into it during an endothermal process, in order to maintain 
the constant equilibrium temperature. 

In the case of a periodical stress (8), p = p,sinft = P, sinw6, 

w = pi/p[0-dpw sin(2w6+ «)— (1—7)w(1+w?)— sin w6 e-F+- 

+0-5(1—r)w?(1-+-w?)-1]. (22) 
The power expended is composed of three parts. The first term is the elastic 
part and disappears in a cycle. The second term represents elasticity too, 
as it disappears in a cycle, but decreases exponentially. The third term is 
constant: it is positive for 7 < 1, negative for r > 1, and vanishes only 
for 7 = 1. When positive (7 < 1), it represents viscous dissipation. But 
what is the meaning of a ‘negative’ dissipation (r > 1) ? 

This seems in flagrant contradiction with the second law of thermo- 
dynamics, according to which in closed systems heat may not be extracted 
from the body and transformed into work. This seeming contradiction may 
be removed by one of the following hypotheses: 

(a) The body does not undergo any internal structural changes. In this 
case we have always ss ete (23) 
hence, by (4), a <<, &; G< &. (: 


The time of relaxation R may not exceed the time of retardation L, and the 
rigidity modulus G may not exceed the firmness modulus H. 

(6) It is a matter of common observation that many bodies, when sub- 
mitted to periodical shear stresses, are gradually destroyed and undergo 
internal structural and physico-chemical changes. Assuming the possi- 
bility of such changes, which may be accompanied by the liberation of 
sufficient amounts of energy, a ‘negative’ dissipation may be allowed for. 
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It is further assumed that the rheological coefficients vary but very slowly, 
so that the above formulae maintain approximately their validity. In this 
case T may exceed 1, R may exceed L, and G may exceed H. Such bodies 
may be called endous solids, by analogy with endothermal chemical 
reactions. 

Bodies with 7 < 1 may be called exothermal or dissipative bodies, as they 
lose mechanical energy which is transformed into heat by viscosity. Here 
belong Maxwell’s firmo-viscous fluid (7 0, G 0, L+o), Kelvin’s 
elastico-viscous solid (7 0, H +o, R= 0), and all bodies with 7 < 1. 
Such bodies may undergo structural changes as well, which mask the 
phenomenon. 

Bodies with t = 1 may be called homothermal bodies. The dissipative 
term is zero, yet there may be a viscous dissipation compensated by some 
change in the structure of the body. This is the basic difference between 
these homothermal solids, both viscous, firm, and elastic, with relaxation 
andretardation, yet without any permanent strain; and the ordinary elastic 
Hookean solids which neither show permanent strain nor undergo any 
structural change. 

The most extreme case of endous bodies is that of the hypothetical pure 
endous solid (r > 00, Goo, H +0). Its rheological equation is 

€ np. 
The power expended is 
w pe npp 0-5 6?( p?)/ot?— np. 
The first term disappears in a cycle; the second term is always negative. 
may therefore be termed endosity, as if it were some kind of negative 


viscosity 


7, Classification of linear bodies 

The above considerations lead towards a new classification of the linear 
bodies, according to the value of their time factor r: 

(a) Exothermal or dissipative bodies (r <1). Energy is dissipated by 
viscosity and transformed into heat. They warm up adiabatically. They 
we more firm than rigid (H > G). Lagging is more pronounced than 
relaxation (L > R). The strain is smaller than in Hookean solids of 
same G. One important case makes a subdivision of its own: 

(b) Maxwell fluids (r = 0). They have no elastic rigidity under constant 
stress (G' = 0). They are firmo-viscous, as they show seemingly elastic 
firmness (H) under periodic, impulsive, or transient stresses. Their rheo- 
logical equation is 


a,pt+a,p=—a,é or p+Rp=—pé or p/p+p/H =. 

















408 S. IRMAY 


Here belong the usual Newtonian viscous fluids having infinite firmness, 
Their equation is ; ‘ 
a,p=—aé or p= ype. 

(c) Endous solids (r > 1). Energy is liberated by internal structural 
changes, the final result being that of a ‘negative’ dissipation. These 
bodies cool down adiabatically. They are less firm than rigid (H < @), 
Lagging is less pronounced than relaxation (L < R). The strain is greater 
than in Hookean solids of same G. Owing to structural changes these bodies 
probably disintegrate after numerous cycles of stress (fatigue). 

(d) Homothermal solids (r = 1). Energy is neither dissipated nor ab- 
sorbed in a cycle, as the viscous dissipative effects and the endous effects 
(due to structural changes) compensate each other. An adiabatic process 
is also isothermal. Lagging equals relaxation (L = R). Rigidity equals 
firmness (G = H). The strain is that of the Hookean solid of same G. 
Their equation is 

pt+Rp = Gle+ Re). 


The Hookean elastic solid (R, L, wu, n = 0) is a particular case 


a,p=a,« or p= Ge. 


(e) Ideal bodies: 
(i) Inviscid ideal fluid (u, R, G = 0), 


ap=0 or p=20. 
(ii) Rigid ideal solid (L, n = 0; G>o), 
age = 0 or e== @, 


All the above bodies are summed up in a table. 




















Body | a, | a | as | a | r | Gi 2} & | I | Tn) 
1. Exothermal or dissipative bodies (r < 1) | | } | 
(a) Fluids (G = 0) | | | 
(i) Newtonian (viscous) . . | ! o|o | | o|ola | o | ao] p 
(ii) Maxwellian (firmo-viscous) . | 4 | 4 ° | lo|lol|dH | H co | pw |e 
(b) Kelvin solids (elastico-viscous) . | 4 | o + | 4 | o |G om} 0 | L 1 GL je 
(c) General exothermal solids T } G |G/r | Gn } 7G 7| nG*/r o] 
2. Homothermal solids (r 1) | | | } } 
(a) Hookean elastic solids. . .|/+]o0/+]o}] |G] o | o | o jo 
(b) General homothermal solids. EA + }+]4 1|G | G | Gn Gn | 7G? | 9 
3. Endous solids (+r 1) ; . ; ! +} +] 4 t | G |G/r| Gy | 9G/r| nG/r | 9 
4. Ideal bodies | | 
(a) Inviscid fluid. : ; aie o|lo| ° ° | o 


(b) Rigid solid 
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SUMMARY 


A solution is derived for the elastic buckling load of a symmetrical I-section 
member subjected to an axial thrust applied through the centroid combined with 
unequal terminal bending moments acting in the plane of the web. An energy method 
is employed, using an estimate of the deflected form obtained by a method of successive 
approximation. The solution may be obtained to any required degree of accuracy, 
Approximate formulae which always give conservative results are also derived. 


1. Introduction 

THE behaviour of symmetrical I-section members under centrally applied 
axial loads combined with terminal moments acting about the minor axis 
has been described by Pippard and Baker (1). The elastic buckling load is 
equal to the Euler crippling load for the member treated as a pin-ended 
strut, but the column will actually collapse at a load less than this since the 
terminal moments will give rise to deflexions perpendicular to the web, 
thereby producing at some stage below the Euler load extreme fibre stresses 
greater than the yield stress. Expressions giving the maximum fibre stresses 
are quoted by Pippard and Baker. 

When the terminal bending moments are in the plane of the web a different 
phenomenon arises. The elastic buckling load is reduced below the Euler 
crippling load due to the tendency for failure to occur by bending about the 
minor axis combined with twisting. The mode of failure is similar to that 
of I-section beams bent about the major axis and failing by lateral instability. 
A solution for columns of rectangular section when the terminal moments 
are equal and bend the column in the same sense has been given by Timo- 
shenko (2), and a similar formula has been derived for symmetrical I-section 
members (see Fig. 1 (a)) by Johnston (3). It is assumed that the ends of the 
column are restrained against twisting about the longitudinal axis, but are 
free to rotate about the minor principal axis. It may be noted that the 
problem of a column loaded as above may also be interpreted as that of a 
column with an eccentric end load acting in the plane containing the minor 
principal axis (Fig. 1 (b)). 

The solutions given in references (2) and (3) are approximate and ignore 
terms which may assume considerable importance. In particular, no 


[Quart. Journ. Mech. and Applied Math., Vol. VII, Pt. 4 (1954)] 
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account is taken of the reduction of torsional rigidity due to axial load, 
in effect noticed by Wagner (4). General equations governing the 
behaviour of members of open section subjected to combined axial loads, 
twisting moments, and bending moments have been presented by Goodier 











: P 

| “eb 

i iY 

| | | 

a) (b) 

\ 4AM 4 A 

| +e 

P Pp 


Fic. 1 


5), Timoshenko (6), and Bleich (7), and from these equations a solution 
to the problem shown in Fig. 1 (equal end moments) may be derived. Even 
the solutions so obtained do not allow for the secondary twisting and 
bending moments attributable to curvature about the major axis just prior 
to buckling. These secondary moments influence the elastic buckling load, 
and their effect has been described in the case of members with no axial 
load by Davidson (8). The influence of these moments is small, however, 
and the general validity of the Goodier solution has been confirmed in tests 
carried out by Hall and Clark (9). 

Solutions for columns subjected to unequal end moments M, and M, (see 
Fig. 2) have been given by Salvadori (10). These solutions were obtained 
by an energy method, correct numerical results being given for ratios of 
end moments of 1-0, 0-5, 0, —0-5 and —1-0. The results are given in the 
paper, in both tabular and graphical form, for a range of values of the axial 
load and discrete values of the warping rigidity factor. Results obtained 
by interpolation are also quoted for other ratios of end moments. No 
formula is given from which an accurate numerical result may be computed 
for any given axial load, warping rigidity, and ratio of end moments. 

The present paper expounds in section 3 a method of solution also based 
on energy equations, but differing from that employed by Salvadori. 
A process of successive approximation is first applied in order to obtain 
a close approximation to the deflected form, this method being of a general 
nature so that any number of terms in the Fourier analysis may be taken. 
This process is much accelerated by using as the first approximation to the 
deflected form values obtained from graphs given in the paper. The energy 
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equation then gives the required result, which may be obtained to any 
required degree of accuracy. When only approximate results are required, 
it is possible to use simple formulae which, when in error, always give 
results on the safe side. Such formulae are discussed in section 4. _ 
The present paper is concerned entirely with theoretical elastic buckling 
loads, and it is assumed that the maximum stresses at the critical loads are 
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less than the yield values. The actual collapse loads of members of low to 
moderate slenderness will therefore be less than the values given by the 
derived formulae. It is also assumed in the analysis that the members are 
initially free from imperfections, the presence of which would cause even 
very slender members to collapse a little below the values given. The 
solutions are presented as the necessary basis for further investigations 
rather than for their use in giving directly the collapse loads of real members. 


2. Notation 


E Modulus of elasticity. 

G Shear modulus. 

d Distance between centres of gravity of flanges of an I-section 
member. 

A Area of cross-section. 

hy t, Moments of inertia about XX (major axis, perpendicular to 
web) and YY (minor axis, in plane of web) respectively. 

I, Polar moment of inertia (J, = [,+-1,). 

Vig By Radii of gyration about XX and YY respectively (A A=, 

Ars = 1). 
‘. Polar radius of gyration (Ar? = I,). 
K Torsional constant of I-section member. 


Flexural rigidity about minor axis (B = EJ). 
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0 Torsional rigidity (C = GK). 

y Warping rigidity factor for member (y = 77d?B/4I?C). 

P Centrally applied axial load. 

D Mean axial stress (Ap = P). %: 

’ Factor giving reduction in torsional rigidity due to axial load 
(a pl,/C). 


M,,M, Terminal bending moments about major axis (Fig. 4). 

8 M,/M,. 

R Shear force in plane of web just prior to buckling (Fig. 4). 

T) Torque at ends of column during buckling (Fig. 5). 

M.,M, Bending moments acting about major and minor axes of a 
cross-section respectively during buckling. 


T Torque acting at a cross-section during buckling. 

P, Euler crippling load for column failing about minor axis as a 
pin-ended strut (P, = 7?B/I*). 

M, Critical value of uniform bending moment which, when applied 


about the major axis, would just cause lateral instability, the 
ends being free to rotate about the minor axis, and resistance 
to warping being ignored (Mj, = 7? BC/I*). 

Z Distance along column from origin at one end. 

u Displacement during buckling of centroid of section in plane 
OZX (Fig. 5). 

v Displacement in plane OY Z. 


Y Angle of twist during buckling about axis OZ. 

Up Strain energy due to flexure about minor axis. 

Us Strain energy due to differential flexure of flanges about minor 
axis. 

U, Strain energy due to twisting. 

Up External work done by axial load P during buckling. 

Uy External work done by terminal moments during buckling. 


3. Solution for the general case 


Consider a symmetrical I-section column of length / subjected to an axial 
load P and unequal terminal moments M, and M, applied about the major 
axis (Fig. 2). Let M,//, = 8, so that all possible cases will be included if 8 
is allowed to vary between —1 and +1. Just prior to buckling the con- 
ditions will be as represented in Fig. 4. Take rectilinear axes OX, OY, OZ, 
the origin being at the top of the column, and the web lying in the plane 
0YZ. The longitudinal central axis is OZ, along which the axial load P acts. 
The moment M, is applied about the axis OX, and M, about the parallel 
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axis O'X’. Equilibrium is maintained by shear forces R acting at O and 
in the plane of the web, where 


Rl = M,—M, = —(1—f)M 


1° 


The ends are fixed in position in the planes OXY, O'X’Y’, although they 
are free to approach each other in the direction OZ. The end sections are 


al 
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free to rotate about the minor axes OY and O’Y’, but are restrained from 
twisting about the longitudinal axis. 

When buckling occurs, the column will take up a slightly deflected form 
under neutral equilibrium, as shown in Fig. 5. Let the centroid g of any 
section distant z from the origin O deflect wu in the direction XO and vin 
the direction YO. Let the same section rotate through an angle @ about the 
axis OZ according to the right-hand rule. There will in general have to be 
equal and opposite torques 7, applied at the ends of the column about the 
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and 0 | axis OZ in order to prevent twisting of the end sections (see Fig. 5), 7, being 
regarded as positive at O when acting according to the right-hand rule. 
Let the bending moments about the major and minor principal axes and 


h they y 


MS are xX 








1 from the twisting moment about the longitudinal axis at any section distant z 

irom the origin be denoted by M,, M,, and T' respectively. The values 
d form of M,, M,, and 7 are again taken to be positive when they act according to 
of any the right-hand rule on that part of the member remote from the origin O. 
id v in Expressions for M,, M,, and T may be obtained from equilibrium con- 
ut the ‘iderations. It will be assumed that the moment of inertia J, about the 
» to be major axis is considerably greater than that about the minor axis, J,, and 


ut the lence in establishing the equations of equilibrium it will be assumed that 
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v = 0. Ignoring second-order terms, it is thus found that 


M, = u,{1\—a—p)}}, (1) 
M, = —Pu—M,{1—(1—A)}] 8, (2 
7 z d , 
T = Ty—M,{1—(1—p)3| M(1—p)*. (3) 
—) 





Direction of ——_| 
longitudinal fibres 
after twisting 











Component of 
longitudinal force 
on area dAin 


plane OXY 
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The bending moment MV, causes flexure about the minor principal axis, 
giving the relationship is - du 
M, = tly 
When considering the relationship between the torque 7’ and the angle 
of twist per unit length (d6/dz), it is necessary to consider the reduction of 
torsional rigidity due to axial load and the effect of differential bending of 
the flanges (warping rigidity). With regard to the former, if a member 
sustains a mean longitudinal stress p and undergoes a twist per unit length 
of the longitudinal axis OZ of d0/dz (Fig. 6 (a)), the stress p acting on an 
elemental area 5A of the cross-section will have a tangential component 
of (p 8A )r (d8/dz) (see Fig. 6 (b)). This will cause an integrated torque about 
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(4) 
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the axis OZ for the whole section of > pr? (d0/dz)5A = pl, (d0/dz), where I, 
is the polar moment of inertia (J, = [,+T,). Hence the effective torsional 
rigidity is reduced by the axial load from GK to (GK—plI,), where G is 
the shear modulus and K the ordinary torsional constant for the section. 
This effect has been discussed by Wagner (4). In an unsymmetrical section 
the torsional rigidity is also affected by bending moments, but bending 
moments have no such influence in symmetrical sections. 

The effect of warping rigidity may be stated as follows. The contribution 
of the web to the moment of inertia about the minor axis may be ignored, 
whence each flange has a moment of inertia about the YY axis of $/,. The 
resistance offered by the flanges to differential bending therefore adds a 
resisti:.g torque of —}#I, d*(d*6/dz) (see Timoshenko (2)). The symbol d 
denotes the distance between the centres of gravity of the flanges (Fig. 3). 
The torque 7’ at any section in the column (Fig. 5) is thus related to the twist 
per unit length by the equation 


7 2 3 
. (GK— pl,) « El,d d°6 


dz 4 de ™ 


Consider now the deflexions wu of the column in the OZX plane (Fig. 5). 
Let these be expressed by the half-range Fourier series 


y= > (« »sin"). (6) 


Similarly, let the angle of twist @ be expressed in the form 


“_< ; z 
rf) 1M, S (> sin es ; (7) 
C _~ l 
where C = GK. Now let the torque 7' be eliminated between equations 
3) and (5), and the value of uw from equation (6) be substituted. Then it 
isfound that, if B= El,,a —?,y—7 72 
) C 4 PO 
Mee Mh 8 
dz 7 dz C 6UC 
where 
H = (1— 8) > (a sin") } m1—(1 : -B);| 2, (na, cos Th (9) 


The value of 7, which appears in equation (8), may be derived by 
integrating this equation between the limits z = 0 toz = 1. Hence it is 


for 
ind that T, = 7M, > (na, ,), (10) 
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where 
n odd, : - 


4 (1—B) 
boa ay 


nmeven, x, = ( 


The function H may be expressed as the Fourier series 





D x z 7? 
H — Zz | > a,{(1—B)Anm- T NX} }eos + °, (12) 
m=1'n=1 l M, 
where 
4 
(n+m) odd, A,» =--—<; = 5 | 
a7 (n*—m*) (13) 
(n+-m) even, A,,, = 0 
and eee 
4 _ 2). 2 
(n+m) odd, Xj,,, = —s : ee) 
a2 = (n®—m?)? 
(n+m)even, n~m, A», =O . (14) 
: (1 
n= Mm, A = +8) 


mm — 


The values of @, 7), and H given by equations (7), (10), and (12) respec- 
tively may now be substituted in equation (8). The coefficients of cos(mzz/I) 
may then be equated, and, using the relations (11), (13), and (14), it is 
found that 
l 8 ao = _ n3 1° th, | (15) 
{(1—a)- 192 2} m Lu \(n?—m? ef 2 " 


n 


Im lest 











where the summation ’ includes only those values of n for which (n-+-m) 


n 
is odd. 

The coefficients 5,, in the Fourier series for the twist of the column 
(equation (7)) may thus be derived from assumed values of the coefficients 
a,, in the Fourier series for the deflected form (equation (6)). The coefficients 
b,, may in turn be used to derive an improved estimate of the deflected 
form as follows. 

The minor axis bending moment J, is eliminated between equations 


(2) and (4), and the value of @ given by equation (7) is also substituted. 


Hence ee . 2 
d um P E a eS N. (16) 
ae PF, 1 \ Mp, 
where < / diaand 
N = — (1—B) iI y (On sin **). (17) 
In equation (16), 7 Bil? and M, = 7v(BC)/l. The symbol Py 


denotes the Euler nisin load for the column failing about the minor 
axis as a pin-ended strut. The symbol M,, denotes that uniform bending 
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moment which, when applied about the major axis in the absence of any 
axial load, would just cause elastic buckling due to lateral instability, 
resistance to warping being neglected. 


The function NV may be expressed as the Fourier series 


N > > bn nm |sin we, (18) 
n=1*m=1 
where 8 (1—B)nm 


(7 m) odd, Lam >= ss ST 
ar (n*—m*)* 


(n+m)even, n~AM, py», = O >. (19) 


n m, pie 


bo 


It is now assumed that 


“_\ . NZ 
u sr) 2 la, sin” 
n=1 ~ 


20 
My ; (20) 





Substituting this in equation (16) and equating coefficients of the general 
term sin(n7z/l), it is found that 


a, 2 : | = (1 B)n > —_ = © | +f b,\, (21) 
(7° P P, \iar- — (n* m*)* 2 | 


where the summation 5’ includes only those values of m for which (n-+-m) 





is odd. The improved coefficients a}, may now be used in place of the 
coefficients a,, and the process repeated until a sufficiently close approxima- 
tion to the deflected form has been obtained. 

The deflected form of the column having been derived, the value of the 
critical terminal moment , for a given axial load P may be obtained 
from the energy equation. Let U, and U,, denote the external work done 
by the axial load and the terminal bending moments respectively. Let Up, 
U,, and U;, denote respectively the strain energy due to curvature about 


. the minor axis, differential bending of the flanges, and twisting of the 


section. Then 


l 
. Pf (du\2 
Up == | (7) dz, (22) 
~ 9 im 
whence, from equation 20), 

a*(M, \? ieee se , 

r= SY mg ee a" 
4 ]2 
Also U,, M)1—(i—p) 700 az, 24 

l 1) ( B)> dz? ( ) 
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whence, from equations (7), (17), (18), and (20), 


Uy = a i)" P, st D. [mean — = )| (25) 








l 
, Bf (d*u\? 
Also Up = | (c3) dz, (26) 
2 J z 
whence, from equation (20), 
: a? M, \* . 
Up = r Gr) Pyl p [ n4a;?]. (27) 
Bd? 176\? 
Moreover.  — — | =. ie 9 
oreovel B . (3 dz, (28) 
0 
whence, from equation (7), 
ss 2/7) 
U,; = a (ar y yPy _ 4p? |. (29) 
Finally, Ur = 8 (eT dz, (30) 
, 2 : dz 
0 
whence, from equation (7), 
: a? /M, \? : 
The energy equation is 
Up+Uy = Upt+Upt Up. (32) 
Hence it is found that 
M, ‘<< uh 2h) ' —_ at 
ae | = §(1] —ax)-+- m2yn2b? 33) 
(32) 2, c a; (n pm) > [i x) +m? y}m?b?, |. ( 





This equation gives the larger terminal bending moment 1, which, in 
combination with the other terminal moment BM, and the axial load P, 
will just cause elastic buckling of the column. The quantity y is the warping 
rigidity factor, while the quantity « introduces the Wagner effect. 

The process of obtaining a solution for a given axial load P and ratio 8 
of end moments is as follows. The coefficients a, in the Fourier analysis of 
the deflected form are first assumed. The coefficients 5,, are then computed 
from equation (15). To facilitate this process the v alues of b,,[ (1—a)+m*y| 
may be obtained by summing the values of a, multiplied by the 
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coefficients in Table 1. Thus, if 8 = 0, and it is first assumed that 
— —e — —— * = 
a, = 1, a, = }, a3 = a, = etc. = 0, then 


b,[(1— x)-+-y] = 0-5000+-0-7205 x 5 = 0-5720, 


b,| (1—a)-+4y] = 0-0450+-0-5000 x 4, = 0-0950, ete. 
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An improved estimate of the deflected form is then obtained from the b,, 
values by means of equation (21). The coefficients by which the b,,, values 
must be multiplied and summed to obtain the values of a,,(n?— P|P,) 
are given in Table 2. The values of b,, and aj, so derived are then substi- 
tuted in equation (33) to obtain the critical value of the larger terminal 
moment M,. Ifa very accurate value of M, is required, it may be necessary 
to repeat the process, using the derived coefficients a’, in place of the 
original assumed values a,. 

A guide to the deflected form of the column for all values of B is given 
in Fig. 7. The curves give the correct values of a,, 20s, 3a3, and 4a, when 
P = Oand y = 0. It is only the relative values of the coefficients a,, which 
are of importance, and the absolute values in Fig. 7 have for convenience 
been chosen so that a?+4a3+ 9a2+-16a2 = 1. The relative values of the 
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coefficients when P + 0 may be obtained very closely by multiplying the 
given values of a,, a3, and a, by (1—P/P,), the coefficient a, being left 
unchanged. It is found that the shape of the deflexion curve is not par- 
ticularly sensitive to the values of the warping rigidity factor y and the 
Wagner coefficient «, and hence Fig. 7 may be used to derive the approxi- 
mate deflected form of any column subjected to any ratio of terminal 
moments and any axial load. It will be found that, by using values of a,, 
derived from Fig. 7, a single application of the process described above 
will give a very close estimate of the critical terminal moments except 
sometimes in the region of 8B = —1. 

The behaviour of a column for which 8 = —1 is complicated by the fact 
that there are then two distinct modes of failure which do not interact. 
Tables 1 and 2 show that if 8 = —1 anda, = 0 for all even values of n, 
then 6, 0 for all odd values of m, and vice versa. There are thus two 
distinct modes, one of which involves only odd deflexion coefficients a,, 
while the other involves only even coefficients. When 

y=0 and P/P, < 0-77, 
the even coefficients lead to the lower critical load, while the odd coefficients 
give the lower critical load when P/P,, > 0-77. When f is in the region of 
-1 but not equal to it, the rate of convergence of the process of successive 
approximations to the deflected form may be slow, but it is found that the 
rate at which successive estimates of the critical moments converge is 
usually satisfactory. 

When £ +1 (i.e. when the applied bending moment about the major 
axis is uniform), the only coefficients a, and b,, which are not zero are the 
first terms a, and b,. Hence it follows from equations (15), (21), and 
(33) that 


A 2 2 
| 4 ig a i (34) 
My} (l—a«)+y} Pr 
When P = 0 and y = 0 (warping rigidity zero), 
M, \? : 
= F, 35 
br, an 


where F is a quantity depending only on f. Similarly, when P = 0 and 
© (corresponding to C = 0 or zero torsional rigidity), 

(37) - F'y, (36) 

where F’ again depends only on f. Values of F and F’ are given in Table 3. 

The dependence of M,/M, on P/P, is illustrated in Fig. 8, which refers 

to cases for which « = 0 and y = 0. The value of P/P, is plotted horizon- 

tally and that of (1/F')(1,/M,)* vertically. Results are given for 8 = 1-0, 
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TABLE 3 

Values of F and F’ 

B | F F’ 
ie) 1*000 | 1°000 
o'9 I°107 1°108 
o8 1*230 1232 
o'7 I°371 1°377 
0-6 1534 |) 1°547 
o°5 I'720 1°746 
o"4 1°964 1-981 
o"3 2151 2°260 
o'2 | 2°4061 2°590 
orl 2°777 2982 
° 3°134 3°445 
—_—@®°z 3°529 3°988 
—o'2 3°806 4°616 
—3 4°428 5°324 
—o'4 4°914 6-092 
“83 5436 | 6874 
—o'6 5°948 | 7°590 
—o7 | 6-433 S121 
—o8 6°796 8-331 
—o'”9 6°905 8-12 
—1'o 6°529 7°52 





0-5, 0, —0-5, —0-95, and —1-0. Equation (34) gives a straight line for 
B = 1-0. When 8 = —1-0, two distinct curves corresponding to the two 
distinct modes of failure are obtained. 


4. Approximate formulae 
The determination of safe approximate solutions to the instability prob- 
lem for straight I-beams subjected to complex loads has been discussed by 
Worthington (11). Equation (34) in the paper by Worthington shows that, 
if M,, is the critical value of the larger terminal moment J, for a given 
value of 8B when P = 0, then a safe value of M,, is given by 
Mj, = (F+F'y) Mi. (37) 
Using this equation together with equation (27) in the paper by Worthing- 
ton, a safe estimate for critical combmations of the axial load P and 
terminal moments /, and BM, is given by 
_U+y)_ Gr) i ia ay (38) 
(F+F’y)\Mg} (l—a)+y} °° Pr 
When 8 = 1-0, F = F’ = 1-0, and equation (38) reduces to the rigorous 
solution represented by equation (34). When f is positive, equation (38) 
gives results which are quite close to the correct solutions. The error 


increases as 8 approaches the value —1-0, but results are always on the 
safe side. 
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Fie, 8. Values of 5 M,) when y = 0, a = 0. 


Asimpler formula, which is, however, more conservative, may be derived 
as follows. The value of «/y is 4l?pJ,,/7*d*.B, which may also be expressed as 
P/P,)(2r,,/d), where r, is the polar radius of gyration. Geometrical con- 
siderations show that for all I-sections r, will be approximately equal to 
or less than 4d. Thus for a 12 in. x5 in. x 32 lb. British Standard Rolled 
Steel Beam, r, = 4-94 in., while $d = 5-72 in. Since P/P, < 1, it will 
therefore be safe to assume that a/y < 1. Since, moreover, F’ > F in equa- 
tion (38) (see Table 3), it follows that safe values for critical combinations 
of P and M, will be given by 

(MOP l FL (39) 
\My) FP Py 
This equation differs from (38) in that it entirely neglects the increase of 
carrying capacity due to resistance to warping and the decrease due to the 
Wagner effect. It gives a conservative answer since the former always 
more than compensates for the latter. 
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5. Effect of flexure about the major axis 

In the above treatment no account was taken of the effect of curva- 
ture of the column about its major axis. Deflexions in the OYZ plane 
(Fig. 4) have two contrary effects. The curvature about the major axis 
increases the stability when the axial load is small, the increase in critical 
moment under zero axial load being in the ratio of 1 :[1—(J,/I,)]-* when 
B = +1 (see reference (8)). When the axial load is high, this effect js 
counteracted by the increase Pv in the bending moments about the major 
axis, where v denotes the deflexion in the direction YO. Both of these 
effects will be small in beams with a low ratio of [,/I,. Beams in which 
I, /I, is higher than about 0-2 do not usually fail by elastic instability under 
combined axial load and bending moment about the major axis, since the 
critical moments are found to produce stresses in excess of the yield value, 
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WATER WAVES GENERATED BY OSCILLATING 
BODIES 
By F. URSELL (Trinity College, Cambridge) 
Received 5 November 1953] 


SUMMARY 

A two-dimensional wave motion is set up in an ideal fluid under gravity by an 
infinitely long cylinder with horizontal generators which is oscillating vertically with 
constant frequency and small amplitude. Surface tension is neglected. At a distance 
from the cylinder the motion approximates to a wave train travelling away from 
the cylinder. The amplitude at infinity depends on the shape of the cylinder and 
on the frequency, and this dependence is relevant to the design of wave-makers 
and to the wave-damping of ships’ motions. Approximate methods for calculating 
it have been given by several writers, but their results do not agree with a recent 
exact calculation for short waves. A critical comparison between exact and approxi- 
mate methods is now made for very long and very short waves. The rigorous 
method for short waves is very laborious, and is replaced here by plausible assump- 
tions about the velocity potential. The principal result is that short waves generated 
on one side of the body cannot pass under the body and therefore cannot contribute 
to the wave amplitude on the other side. This is intuitively obvious but is not ob- 
tained by the earlier methods. Green’s theorem on potential functions plays a 
central part and is used to obtain the leading term in the wave amplitude when the 
wavelength is short. The representation of bodies by a source distribution is 
criticized. 


1. Introduction 

WHEN a body floating on the surface of a fluid under gravity is given a 
forced vertical motion (heave) of constant amplitude and frequency, it sets 
up a wave motion in the fluid. At any point of the fluid a stationary state 
is ultimately attained which depends on the shape and size of the body, the 
amplitude of oscillation of the body, and on the frequency. The calculation 
of this final state may be reduced to a problem in classical hydrodynamics 
by simplifying assumptions which are satisfied in many applications. 
Quadratic terms in the equations of motion are neglected; this is justified 
if the amplitude of heave is small compared with the dimensions of the 
body and with the wavelength. The viscosity is neglected since the motion 
atany point is oscillatory, and (ina first approximation for small amplitudes) 
the vorticity is thus likely to be confined to oscillatory boundary layers of 
negligible thickness near solid and free boundaries. (A velocity potential 
therefore exists.) Surface tension is neglected; this is justified if the wave- 
length is large compared with the critical wavelength (1-7 cm. for water). 
Finally, it is assumed in the present paper that the motion istwo-dimensional. 
This will be the case when the wave-making body is an infinitely long 


(Quart. Journ. Mech. ana Applied Math., Vol. VII, Pt. 4 (1954)] 











428 F. URSELL 
horizontal cylinder or a horizontal cylinder confined between walls normal 
to its axis; the velocity component parallel to the axis then vanishes. Many 
of the results obtained here can be extended to apply to three-dimensional 
motions. 

When the foregoing assumptions are made, the problem of determining 
the motion becomes a boundary-value problem of linear potential theory, 
of which the equations are given in the next section. The waves at a distance 
from the cylinder are then regular wave trains spreading out from the 
cylinder, with an amplitude proportional to the heaving amplitude. The 
ratio of wave amplitude to heave amplitude (called the amplitude ratio in 
the present paper) is a function only of the ratio of body dimensions to 
the wavelength, and is important in determining the damping of the free 
heaving motion (1, equation 34). Much of the present paper will be con- 
cerned with its calculation. 

Mathematical techniques for calculating the wave motion (and in par- 
ticular the amplitude ratio) have been described in earlier papers (2, 3). 
For each frequency an infinite system of linear equations in an infinite 
number of unknowns must be solved; the amplitude ratio is a function of 
the unknowns. For slow motions this method is satisfactory, but for high 
frequencies the system of equations is ill-conditioned, i.e. a large number 
of equations has to be solved to a high degree of accuracy in order to obtain 
an adequate approximation to the amplitude ratio. So far computations 
are available only for the half-immersed circular cylinder for values of the 
ratio diameter : wavelength from 0 to 1-5. The amplitude ratio rises from 0 
toa maximum value 0-88 at 0-6, and then decreases to 0-67 at 1-5 (2, Table 1). 
For higher frequencies only approximate theories (e.g. that of Holstein, 4) 
have been available until recently, according to which the amplitude ratio 
has an infinite number of zeros or of small minima as the frequency increases. 
The first minimum occurs in this theory when wavelength and diameter 
are comparable, but a recent exact mathematical investigation on the semi- 
circle has shown that for large frequencies the predicted minima do not 
occur (5, equation 5.7). A critical discussion has thus become necessary 
and is given in the present paper. Since a rigorous treatment involves much 
mathematical detail (as in 5), an alternative course is pursued here: 
plausible assumptions are introduced in the place of mathematical proofs, 
whereby it is hoped to make the discussion available to a wider class of 
readers. 


2. The equations of motion 
The origin of rectangular Cartesian coordinates is taken in the mean free 
surface (y = 0, where y increases with depth). If the cylinder has a vertical 
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WATER WAVES GENERATED BY OSCILLATING BODIES 429 
plane of symmetry, this is given the equation x = 0; otherwise the origin 
of xis chosen arbitrarily. The beam (= length of z-axis intercepted by the 
cylinder) is denoted by 2a. Polar coordinates are defined by the equations 

x = rsin86, y = reos@. 
The time-dependence of the velocity potential $*(x, y, t) is simple harmonic, 
with period 27/o, and it is convenient to use a complex-valued potential 
function 4(x, y)e-*® of which the actual velocity potential 4*(z, y,t) is the 
real part. It is easily seen that 
d(x, yje* = H*(x, y, t)+16*(zx, y, t+-77/20). 

The argument of the present paper is greatly simplified by the use of 
time-independent functions like ¢(x,y). The depth of the fluid is assumed 
infinite; the modification for finite depth is straightforward. 

The equation of continuity 

Ah ad 0 
a2 ' Oy® 
issatisfied throughout the fluid. The boundary conditions are that on y = 0 
; C 
o%p+g% = 0; (2.1) 
cy 
that on the curve occupied by the cylinder in its mean position 
Cc d . r 
—ila cos§ = U, cos, (2.2) 
on 
where 5(2, y) is the angle between the vertical and the normal to the curve 
at any point, / is the amplitude of heave, and @/én denotes normal differen- 
tiation into the fluid; and that at infinity the waves spread outwards: 


d(x, y) > A,e-Kv+iKizi as x—> +00, 


d(x, y) > A_e~Ku+tKit! ag xz—> —oo, 
where 4,, A_ are complex constants (A, = A_ for symmetrical bodies) 
and K = o*/g. For a derivation of these equations see (2). An auxiliary 


potential function needed later represents a wave source at (&, 7): 


G(x, y; €, n) G(E,; 2, y) 





a Oo 2 « OSs k =— 
(x—£P?+(y—n)? ny +m C08 K(Z—E) gy 
): k—K 
— 2mie— K+” cos K (x—&) 
see 5, equation 3.2), which satisfies the equation of continuity, the free- 
surface condition (2.1), and also the radiation condition since 


> —Darie—Ku+nei K\x-€ 
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F. URSELL 
3. Green’s theorem 

The argument is based on Green’s theorem (6, §§ 44, 60): 

If two harmonic functions F,(x,y), F(x, y) are regular inside and on a 
closed curve, then the line integral along the closed curve 

ay ay 
$ EG y) fale, ¥) — F,(z, ne) ds(x,y) = 0, 
on on 
where @/én denotes differentiation along the normal to the curve, and ds(x, y) 
is a line element of the curve. 

Let this theorem be applied to the potential 4(x, y) of the motion, and 
to the source potential G(x, y; €,) where (€, 7) is inside the fluid, and let 
the line integral be taken along the boundary of the whole region occupied 
by the fluid. The neighbourhood of (€, 7) must be excluded by a small circle 
which makes a finite contribution to the integral (see 6, § 57). There is no 
contribution to the integral from y = 0 where both ¢ and G satisfy (2.1); or 
from infinity where both ¢ and G are outward-travelling waves. Thus 


, ; ; 
$(€,9) = — Gee, y; &, 0) (x, y) bla, y) Glo, y; &, n) | ds(a, 9), 
27 J | on on 


where this integration is along the boundary of the cylinder only and the 
differentiation is into the fluid. In words, the motion at any point may be 
thought of as due to a source distribution on the cylinder of strength 
= =f (a, y), 
together with a distribution of normal dipoles of strength —(1/2z7)d(z, y). 
The potential at a distance from the cylinder is found by substituting for 
G(x, y; €, n) its limiting form (2.3), whence, as € > +00, 
d(é, n)> — je— Kn+i KE | [ ~Ky- iK2o (x, y)— A(x, y) ah (e-Ky ixs)| ds, (3.1) 
; on én 

so the wave amplitude is proportional to the last integral. To calculate 
this we must know both é¢/én (the source strength) and ¢ (the dipole 
strength) on the cylinder. The first of these is given by equation (2.2), and 
the second is thus determinate in theory. Our problem is the practical one 
of actually finding ¢ sufficiently accurately to give a satisfactory approxi- 
mation of the integral in equation (3.1). Previous discussions have not been 
based on Green’s theorem but on intuitive physical considerations. The 
most detailed is that of Holstein (4), who from a somewhat obscure argu- 
ment (inspired by Prandtl) concludes in effect (4, equation 6) that the 
dipole term in (3.1) may be neglected.+ He states, correctly, that this 
approximation is valid when Ka is small. 

+ A representation of the motion by sources alone is possible (12), but then the source 


strength depends critically on the frequency for large Ka and can only be found by solving 
an ill-conditioned integral equation. 
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WATER 


4, The semicircular cylinder 


On Holstein’s theory the amplitude ratio of a heaving semicircular 


evlinder (not given by Holstein) is, from (3.1) without the dipole term, 


2Ka | cos @cos(Kasin 0)e~ *4 ©8 & dg 


0 


Ka | Ka . 
; ; -> 9 | Ssinu or sin u 
sin Ka+ Kacos Ka+ K2a? | —— du = —2K*%a* | —— du 
3 
! u J ou 
x x 
2cos Ka l —— 
ol - 3] when Ka is large. 
Ka K*a? “ 


The equivalence of these expressions may be shown by expanding in 
powers of Ka.) This function vanishes for a value of Ka between }z and z, 
and thereafter its absolute value has a minimum (zero) when Ka is an odd 
multiple of $7, and a maximum when Kaisan even multiple, approximately. 
These oscillations in the amplitude occur because the waves are in effect 
generated near the two highest points of the immersed body, and the 
resultant amplitude is due to their interference. Holstein discusses at length 
acylinder of rectangular section (see section 5 below), where interference 
ilso occurs, and while observing that his argument is valid only for small 
Ka he clearly considers the interference effect to be real (4, top of p. 109). 

There is, however, a difficulty here which is striking when the wavelength 
isvery small. The interference between the two sides of the cylinder might 
then be expected to be small, as the only communication between the 
sides is under the cylinder, at a depth of many wavelengths. From this 
consideration the amplitude might be expected to be non-oscillatory for 
large Ka. If this is so, then the dipole term (so far neglected) is as important 
as the source term in determining the amplitude for large Ka. To settle 
this point the case of the semicircular cylinder was worked out, and the 
amplitude was actually shown to be non-oscillatory for large Ka (5, 


equation 5.7). The leading term in the amplitude ratio for large Ka will 
now be obtained by a simpler method using Green’s theorem and plausible 
assumptions, although some parts of the calculation will appear artificial 
when separated from the complete treatment given in (5). 


It should be kept in mind that for short waves the free-surface condition 


2.1) tends formally tod = 0. As Ka + o we therefore expect the potential 
tote.d to the limit ® ,a*r-} cos 6 (satisfying this boundary condition 


and (2.2)) except in a narrow surface layer of thickness 2K-! where there 
re waves. It is known (2, p. 223) that the potential 


. ,({cos@ . cos2é@ 
l gt “TT T}9 
r K ae 
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satisfies the free-surface condition exactly; it also satisfies the condition 
(2.2) on the semicircle approximately when Ka is large. Write 


cos we cos 20 1 
= —U,a{—- i 
? ( r | Kr ) Ka" (4) 
Then on the semicircle 
1 _ _ 91, cos 28, (4. 
er 


which does not vanish near A (see Fig. 1) where it takes the value 2U, 

Thus ¢, approximates to a vertical wave-maker near 

‘7 eat A which communicates to the adjacent fluid a hori- 

y zontal velocity 2U,e-**. What happens below the 

B_ surface layer is irrelevant because of the exponential 

decrease due to depth (see equation (3.1)). If, then, 

we apply Havelock’s wave-maker theory (7, see (5.2) 

below) we may expect and do actually obtain the 

right result (see 5, p. 103). But the argument may 

be made more convincing by the use of Green’s 
theorem. 

Denote by ®, the potential satisfying (4.2) and the limit condition ®, = 

on @= +47. Then 


ry Bs new __])\n41 2n+1 
0, = 8U,0 > = 7 cos(2n-+-1)@, 
" (2n—1)(2n+-3)\r 


a 
Fig. 1. 





which is bounded on r = a, and it is assumed that ¢, tends to ®, boundedly 
on the circle as Ka > 0. Apply Green’s theorem to the functions ¢,(z,¥) 
and G(x, y;, 7), and let the integration be along the boundary of the fluid 
lying to the left of ABC. When é - +00, we obtain as in (3.1) 


ost ek iK(é—-a) 
ig, (€ _ 
“hg 
f $,(8 (e -Kr cos 6,iK(a-r sin 9) __ p—Kr cos OpiK(a-r sin Hor a d@+ 
lay or a 
0 


+ i fou 2; y)z “(¢ -Ky-iK(e—a)) _ ¢—Ky-iK(z-a) 


é : ‘ 

‘#1 (a, »| dy. (4.3) 

ox x=0 

In the last ee e~- Kv never exceeds e~ *@, and so this integral is negligible 

if d, and é¢,/éx remain bounded on x = 0 as Ka > oo. So only (4.3), whichis 
im 


| |- K cos ( e —Ka cos §,iKa(1—sin (0) _. 


—ikK sin @ e—Ka 00s 9giKa(1-sin Dd, (8) __ e—Ka cos 6iKa(1-sin 8) | a dé, (4.4) 
or 
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survives, in which @4,/ér —2U, cos 26, and nothing is known about ¢, 
except that it approaches ®,. The exponential e~*« 5? shows that almost 
the whole amplitude comes from the neighbourhood of 6 = 37. The second 
term of (4.4) is apparently an order of magnitude larger than the first and 
third, but, as ¢, is ultimately small near the surface while @¢,/0r remains 
finite, this appearance is deceptive. Another application of Green’s theorem 
is needed. 

Apply Green’s theorem to the boundary of the region to the left of A BC, 
and to the harmonic functions ¢, and e‘*@-@-Ky, Note that the latter 
represents outward-travelling waves at infinity. The contributions to the 
line integral from @ = 47 and from infinity vanish. The contribution from 
BC is exponentially small. Thus 


14,(8) (¢ Kr cos 6piK(r sin 6—a ) e—Kr cos 6,iK(r sin @- a1 a dé 
| - Or \r=a 


> | | K cos 0 e~Ka 00s ¢-iKa(t—sin g (9) + 
| 


LiK sin § e—Ka cos 6p—-iKa(1—sin 6), (@) —e—Ka cos 6p iKatt—sin ) OP a dé (4.5) 
. or 
is exponentially small. The second term is nearly the same expression as 
n (4.4), except for the negative sign. Add the integral (4.5) to (4.4): 
19, (€, n)é iK( a)oK 
> | 2K cos 6 eK © 903 Ka(1—sin 6)¢,(0)+ 
0 


2K sin 6 e~Ke 008 gin Ka(1 —sin 0)¢,(0)— 


or 


2e—Ka c08 (eos Ka(1 sin) Sta do 


an exponentially small term. (4.6) 


In this expression the unknown function ¢,(@) is multiplied by factors 
which vanish near 6 = 47, and if ¢,(@) is assumed to be uniformly small 
near A, the terms involving ¢,(@) are small compared with 

$7 
( Ka cos 6 e~*2 0s ® cos Ka(1—sin 6) dé Uja 
0 


and | Kasin 0 e-*4 ©8 sin Ka(1—sin 6) dé U,a, 
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Ff 
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i.e. they are small compared with (Ka)-!U,a. But the last term in (4.6) jg 


har 
4U,a e~Ka cos P93 Ka(1—sin 0)cos 26 dé 
0 





1 
— —4U,a|— smaller terms 


and is therefore the dominant term. Thus finally, as é > +« 


— — 1 
$,(€,n) > —4tU a e-B1tikE | z+ amaller terms| 
a 


, »(cos@  cos20 ] 
and ee ee CT cc oh erent ohn ee 
¢ " ( pr rad . Ka™ 
The first term vanishes at infinity and so, when £ > +00, 
- boas 1 4 
$(€, n) > —4iU, ae-*0 “¢-0| = tsmaller terms]. (4.7) 
This is the answer obtained in (5). The velocity U, = —ilo, where / is the 


amplitude of heaving, the wave amplitude at infinity is |icg-14|, whence 
the amplitude ratio behaves asymptotically as 


4/Ka, (4.8) 
which is non-oscillatory. The asymptotic behaviour of the phase is also 
determined by (4.7). In the notation of (2), as Ka > o, 

A(Ka)+iB(Ka) = 4rK?a*e"Ke-1*) 4 smaller terms, 


a » ‘ 
whence tan" — Ka+47—>0 as Ka > oo. 


This last result is in particularly good agreement with the computations 
because the expression on the left happens to vanish when Ka = 0. The 
virtual mass can also be found with little trouble. The dynamic pressure 
is ipod, and its resultant force on unit width of the cylinder is 


ia is 7 s 26 ] 
2iapo | $(8)cos 6 dé = —2iapo cos 6| —Uya{ cos 0 + )+z,%| dé, 
0 0 <7 


where p is the density of the fluid. If we replace ¢, by ®,, we have 


Ka 


}cr har 
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snd the virtual-mass coefficientt is 1 —(4/37Ka)+smaller terms. On review- 
ing the preceding calculation we note the essential role played by the 
subsiciary function ¢, defined by (4.1). This function resembles the poten- 
tial due to a vertical wave-maker in the surface when Ka is large, and this 
fact guides the subsequent calculation involving Green’s theorem. In the 
rigorous treatment (5) no use was made of ¢, but a second approximation 
to ¢ was found by a direct method. 

It has been shown that for sufficiently large Ka the wave amplitude due 
toa semicircle is not an oscillatory function of Ka. It seems probable 
that the amplitude ratio (tabulated as Ka(A?+ B?)-* in (2)) has only a 
single maximum 0-88 at Ka 0-6 and then decreases steadily as Ka 
nereases to infinity. This supposition fits in well with the asymptotic 
value 4/Ka and with the computed values. Also, Holstein’s experiments 
4. Figs. 4-7) show no subsidiary maxima for an immersed rectangle unless 
the depth of immersion is a small fraction (less than }) of the beam, whereas 
for a semicircle the effective depth of immersion is more nearly half the 


beam for short waves. 


5. The rectangular cylinder 
This is the case studied theoretically and experimentally by Holstein (4). 
The rectangular cylinder, of beam 2a, is immersed to a depth h (Fig. 2). 


Neglecting the waves due to the dipole distribution, 
~=<-2a 







Holstein finds an amplitude ratio ——-4 


2e-Kh sin Ka, (5.1) 


h 


where K = o?/g (4, equation 11), which has an infinite l 





number of zeros. It will now be shown that the exact B 
theory gives a different answer. Let ABC be the 
vertical through A. c 

Fic. 2. 


By Green’s theorem applied to (a, y) and G(z, y; €, ) 
and integration along the boundary of the fluid to the left of A BC, one finds 
that, as € > +00, 

id(E, n)eXn-*KE-a) 


| | A(x. y) : (¢ Ky—ikK(2£-a ) __e—Ky-iK(x—-a) cr) dy; 
zZ=a 


Cx Cx 


« 


0 


ind, by Green’s theorem applied to ¢ and e~Ky+iK@-a), 





— C 
d(x, y) © (e-Ku+ike a)) _ ¢—Ky+tK(z—a) ceu)| dy = 0. 


+ This coefficient is 8/7? times the coefficient m(Ka) defined in (2). The left-hand side 
f the last equation on p. 224 of (2) should read 2m(Ka)/7. 
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By addition, as € > +00, 


—id(€, ») >» —%e-Kn+ iK(€-a) [ em] 


0 


A en)| dy, (52 
OX @r=a 
a result due to Havelock (7, equation 9). The integrand of (5.2) vanishes 
on AB, the vertical face of the cylinder. If Ka is large, then @4/éx may be 
replaced on BC by é®/éx, where ® is the limit potential satisfying 0 owt 
on on 
on the cylinder, and ® = 0 on y = 0. It thus appears from (5.2) that, in 
this case also, the amplitude at high frequencies is not an oscillatory 
function of the frequency. The amplitude is exponentially small; its exact 
form depends critically on the shape of the cylinder near B. Holstein’s 
experiments do not agree with his formula (5.1) and he finds no minima in 
the amplitude except when h is small 


(h=ja=1-5cem., and h=}a=1cm.). 
The experiments must have been difficult for such small depths, but it may 
be that these minima are real. If this is so, then the flat plate in the surface 
will repay further study. The asymptotic formula may perhaps be found 
by a sufficiently ingenious application of Green’s theorem. 

Holstein uses energy considerations instead of Green’s theorem to find 
the amplitude at infinity. This is a satisfactory method when ¢ is known 
exactly on the cylinder. But actually ¢ is known only approximately, e.g. 
for the semicircle by equation (4.1) with ®, replacing ¢,, which shows that 
to this approximation the dynamic pressure ipod is in quadrature with the 
normal velocity. So to this order no work is done and we obtain merely 
the trivial result that the amplitude is small for large Ka. But by the use 
of Green’s theorem the leading term (4.8) in the amplitude can be obtained, 
Thus the energy method is comparatively inefficient for large Ka. 

Among others who have written on this problem may be mentioned 
Weinblum and St. Denis (1) and Dimpker (8). The assumptions of the 
former appear to be similar to Holstein’s. A different approach is used by 
Dimpker. To solve the problem of the flat plate in the surface he replaces 
the plate by an elliptic alternating pressure distribution acting on the free 
surface. For long waves this leads to the amplitude ratio }7K?a? instead 
of the correct value 2Ka. Nevertheless, an experiment showed that the 
former fitted the observations better than the latter (8, Fig. 5) in the range 
0:8 < Ka < 1-5 when h = ja, see Fig. 2. This shows that the exact 
asymptotic formula 2Ka is of limited applicability, and this is confirmed 
by the computations for the semicircle where the accuracy is within 25 per 
cent. only if the wavelength > 30a, i.e. Ka 2. For high frequencies, 
Dimpker’s theory also predicts an infinite numver of maxima and minima. 
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Concluding remarks 

To sum up this discussion, it is physically obvious that waves generated 
at different points of a floating body cannot interfere if the points are sepa- 
rated by the body and the waves are short. All the unsatisfactory theories 
described in the present paper neglect this shielding effect which may, 
however, be allowed for in simple cases by Green’s theorem (sections 4, 5 
above), but so far only a beginning has been made. A shielding effect may 
arise in a wider class of problems. For instance, Michell’s theory of displace- 
ment ships (reviewed by Guilloton, 9) assumes that a ship of small beam 
ean be represented by a distribution of Kelvin wave-sources only. It is 
found that this theory predicts maxima and minima in the wave-resistance 
which are not verified experimentally. Havelock (10) has suggested a 
viscous mechanism to account for the discrepancy, and has discussed the 
relative efficiency of bow and stern in making waves (11). In the light of 
the present work, we may suspect that a part of the discrepancy is due to 
the finite width of the beam and its shielding effect (represented by a dipole 
distribution), without taking account of viscosity. The limiting case of very 
short waves (= very small speed ahead) may be tractable by Green’s 
theorem, or by an extension of the method of (5), although the asymmetry 
fore-and-aft of Kelvin wave-sources (6, section 256) is likely to cause 
complications. 
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A VORTICAL SINGULARITY IN CONICAL FLOW 


By M. HOLT (Armament Research Establishment, Fort Halstead, Kent) 


[Received 9 July 1953] 


SUMMARY 


The singular line introduced by Ferri in the supersonic flow past a yawed cone 
is analysed. It is shown that such a singularity can arise only at points where the 
velocity is directed along the radial line, that at the singularity the velocity is many. 
valued, the vorticity is infinite, but the pressure is single-valued. The singularity 
can be fitted to a general field of conical flow behind a shock. Flow in its neighbour. 
hood can be regarded physically as that due to a pencil of vortex sheets of varying 
strength. 


List of symbols 


r,@,¢ spherical polar coordinates based on the origin of the conical 
field; 


u,v,w velocity components in the r, 6, d directions respectively; 


p pressure; 
p density; 

a velocity of sound; 

U* limiting velocity; 

S entropy; 

€,7,¢ components of vorticity. 


Symbols relating to Cartesian coordinates are defined in the text. 


1. Introduction 

[N the more important papers written so far on the problem of supersonic 
flow past a yawed cone, two distinct approximate methods have been 
employed. The first (1), due to Ferrari, makes use of the hodograph in 
each meridian plane and is intended to be applicable at any angle of 
incidence. The second, initiated by Stone (2) and developed by Kopal (3), 
is essentially a method of perturbation of axially symmetrical flow past 
a cone treated by Taylor and Maccoll (4). Both methods have been 
criticized by Ferri in an important paper concerned with the physical 
properties of flow past a yawed cone (5). Ferri points out that both methods 
involve the assumption that different members of the family of constant- 
entropy surfaces may intersect the solid conical surface at distinct points. 
He shows that this is not valid since the solid surface itself must be a 
surface of constant entropy, and then proves that all constant-entropy 
surfaces, including the solid surface, must intersect along the generator 
lying in the meridian plane on the leeward side (see Fig. 1). 


[Quart. Journ. Mech. and Applied Math., Vol. VII, Pt. 4 (1954)] 
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A VORTICAL 


For the problem of the slightly yawed cone Ferri estimates the correction 
which has to be made to the results already obtained so as to take account 
of hi: , nysical picture. He represents the effect of convergence of constant- 
entropy surfaces towards the singular generator by the introduction of a 
yortical layer around the surface of the cone. The application of this 
correction certainly yields results in closer agreement with experiment. 

In the general case of large angles of yaw a more fundamental change in 
the existing method will be required, which must take account of the com- 
plete properties of the singularity. Accordingly, in this paper, as a first 
step towards the full solution of the problem of the yawed cone, we 
introduce a singularity which will satisfy the physical requirements laid 
down by Ferri and then describe its more important properties. 


2. Mathematical definition of the singularity 

We employ spherical polar coordinates r, 6, ¢ based on the origin of a 
conical field of steady supersonic flow and chosen so that their axis is along 
the singular line. In such a field no fundamental length is involved and the 
velocity vector, pressure, density, and other physical variables are functions 
only of angular variables. The equations of motion are then independent 
of rand may be written in the form 


ou ou 9 9 
v 1 wcosec G —-v*—u* 0, (2.1) 
oo Oc 
ra ov a Cp 
p) lay w cosec 0 — UV w- cot g\ = —_a2Zt > (2.2) 
“| O6 c d | fal 
; COW Cw P Cp P 
py\v w cosec @ | ww—+t-vw cot | = -a* cosec 0 ot ’ (2.3) 
1 o0 CO od 
ov ow Cp Cp 
py! Qu 4.7 cot 6+-cosec 6 - —vWv A -weosec 6, (2.4) 
00 od) fal] C 
2 1 /, )*2 2 32 2 95 
with a t(y—1}{U**—u?—1 wri. (2.5) 


where u, v, w are the velocity components in the r, 0, and ¢ directions 
respectively, a is the velocity of sound, p is the pressure, and U* is the 
limiting value of the velocity attained when p = 0. It is assumed that a 
polytropic equation of state is satisfied. 

We consider a singularity located at the origin in the neighbourhood of 
which the pressure and velocity components are defined by expansions in 
series of the form 

u Ug(P) + Ou, (h)+-6?u(h)+..., 
v Vo(d) + Ov, ($)+67v,(h)+..., 

W = Wo(d)+ Ow, (d)+6?w,($)+..., 
p Poh) +Op,(b)+6?p.(d)+.... 
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If we substitute these expansions in equations (2.1) to (2.4) and equate 
coefficients of powers of 0, we obtain a series of equations relating w,, »,, 
Wo, Pos Ur» V1, Wy, Py; ete. The terms independent of 6 give 


Wty = 0, (2.6) 
Poy(Wo%o—U) = 9, (2.7) 
Po Y(WyWo+Uy Wo) = —Hy—1)(U*?—up—vp—us)po, (2.8) 
Po Y(Vo+U%) = —W Po, (2.9) 


where uy = Ou,/éd, ete. 

The choice in equation (2.6) leads to two distinct types of singularity, 
The first, defined by w, 4 0, uy = 0, can exist-only in isentropic flow. The 
second, defined by w, = 0, uy) 4 0, arises in non-isentropic flow and is 
therefore the type we require in dealing with the yawed-cone problem. 


3. The isentropic singularity 
Suppose that, in equation (2.6), w) 4 0, uy = 0. Then 
Up constant, (3.1) 
and from (2.7), since p, must be non-zero, we obtain 
Uo = We (3.2) 


Solving (2.8), (2.9), and (3.2) for w», wo, and pp» we find 


Vy = (A/n)sin(nd+e), (3.3) 
Wy = Acos(nd+e), (3.4) 
log py = y(1/n?—1)log{cos(nd+)}-+ constant, (3.5) 


where n? = (y—1)/(y+1), A? = n*(U**—u@), and eis an arbitrary constant. 
From the equation of state and equation (2.5) we obtain the relation 
between the entropy, S, and wu, v, w, p, 


S—S = c,ylog{U*®—uw—v?—w*}—c,(y—1)log p+constant. (3.6) 


If we assume that 
S = S(d)-+08,(4)+...5 
and equate terms independent of 6 in (3.6), we obtain a relation between 4, 
Uo, Vp, Wo, Po. If we substitute the expressions (3.3), (3.4), and (3.5) for 
Up, Wo, log po in this relation, we find 
S,(¢)—S = constant. (3.7) 
Since all constant-entropy surfaces meet in the singularity, equation 
(3.7) shows that the entropy is constant throughout the field. Thus this 
singularity exists only in isentropic flow. 
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4, The non-isentropic singularity 


If we now assume that 
Wo v0, Uy FY, (4.1) 


equation (2.7) is satisfied identically and equation (2.8) gives 


Po = constant; (4.2) 


equation (2.9) then reduces to 








iy 
i ee 
Fig. 1. Section of flow past a yawed cone on a sphere r = constant. AA’ is the 
eridian plane, BC_B’D is the conical shock, E is the axis of the solid cone, and F 


eines the undisturbed stream direction. The lines converging on O are in the 
rection of the 


velocity component tangential to the sphere r = constant. Within 
the shock they are also the traces of surfaces of constant entropy and vorticity. 


If we now equate coefficients independent of @ in (3.6), we find 


So(¢)—S = ce, y log(U**—uj)—c,(y— 1 log py+ constant, (4.4) 
So()—So(0) = c, y log{[ U **—ug(¢) |/[U *?—u2(0)]}. (4.5) 


‘hus the entropy is many-valued at the origin, so that (4.1) indeed defines 
‘singularity, such as that illustrated in Fig. 1, in which constant-entropy 
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surfaces may intersect. We therefore examine its properties in greater 


detail. 


When (4.1) is satisfied, we obtain the following relations between 1,, 


V1, Wy; Py: 


W, Uy = 0, 


ne 


4(y—1)(U*2—u2)p, = 0, 
4(y—1)(U*?—u2)p, = 0, 


») On i J x 
Poy(2Up+2r,+u;) = 0. 


(4,6) 
(4.7) 
(4.8) 


(4,9) 


As in Fig. 1, all the surfaces of constant entropy converge on the singu- 
larity, and, provided that their curvatures remain finite, the final values of 
éS8/e@ must vanish. In Fig. 1 the final curvatures of both bounding surfaces, 
namely the plane of symmetry and the surface of the cone, are certainly 
finite, and we assume that those of surfaces within these boundaries are 


also finite. We assume that the same is true for any singularity of this type. 


We therefore obtain the further relation 


(75) = 0, for all ¢. 
26 } 9-0 
From (3.6), 
os _ : ; 
c6 c6' a0 


00 p 


so (4.10) becomes 


9 = 
2yUo% , (y—l)py — 0, for all ¢. 


U— us Po 
From (4.7), 9, == 0; 
hence, from (4.11), Uy 0: 
from (4.6), w, = 0, 
and, from (4.9), Uptv, = 0. 


Equation (4.8) yields no new relation. 


Tae: ‘ ou ov ow\ c,(y—1) ep 
- — 2c, y(U*2—u2—v2—w?) (n 41 1 o- O,ly- pA 


(4.10) 


ay 


(4.11) 


(4.12) 
(4.13) 


(4.14) 


(4.15) 


We now equate coefficients of 6?, 6°, and 6* respectively in equations 


(2.1) to (2.4) and find: 


From 6? terms W> 0. 





(4.16) 
(4.17) 


(4.18) 
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From #* terms, Ps = 9, (4.19) 
Uy Ws = Up(Uy+ 2), (4.20) 

2u,+3v3+u3+4uU, = 0. (4.21) 

From 64 terms, Uy Wy = 3Uy Us, (4.22) 
Do Y(3Ug Vg + Up Ug + Up Ws) k(y—1)4p,(U*?—u5), (4.23) 

3Py Up Ws = Hy—1)(U*?—Uuz)p;, (4.24) 

2u,+4v,+w, 0. (4.25) 


4.21), (4.23), (4.24) are four 
inear first-order equations for w,, ve, We, and p,. Eliminating wu,, vs, w. 
| 2> “3 3 4 5 “a, “3 3 


If we regard wu, as known, equations (4.20), ( 


from these, we obtain the following second-order linear equation for p,, 
pPy— F' (¢)p4+24p, = u2/K(U*?—u), (4.26) 
where F(¢) log(U*2—u2), and K = (y—1)/2poy. 
Equation (4.26) can be solved for p, when wy is given as a function of ¢. 


3 log u, 


The other three coefficients u,, v3, w, are then immediately determined 
without further integration. Since w, is variable, it follows from equation 
4.26) that p, must be variable also. In general, the remaining coefficients 
will be variable and non-zero. 

The terms in @° yield two equations in p,, us, v4, and w, which, combined 
with (4.22) and (4.25), determine these coefficients. The two remaining 
equations introduce new variables w;, u,, and v;; to determine these we 
require two further equations from 6* terms. A similar procedure is followed 
in finding all higher-order terms. Thus all the coefficients in the expansions 
ean be found, to within an arbitrary constant, once w, has been prescribed. 

Since u)(#) is an arbitrary function, it follows from equation (4.5) that 
§,(¢) is also. Thus the singularity can be introduced into a field of conical 
flow behind a shock wave with any required variation in strength; in par- 
ticular it can be employed in the solution of the yawed-cone problem. The 
In a field of flow 


disturbed by a yawed circular cone such points can arise only in the plane 


singularity can exist only at points where vy, = w, = 0. 


ofsymmetry, where they must lie either on the surface itself or, in the case 
of large angles of yaw, displaced from the surface on the leeward side. 


Expansions for wu, v, w, p in the neighbourhood of the singularity have 


f 


the form 


U = Uyt+u,67-+..., 

9 9 3 

a Uy O-+-v, OF +-..., 
w w, +-..., 


Pp = Pot p,-+..., 
where U., v's, Ws, py are related to uy by equations (4.20), (4.21), (4.23), and 


(4.24), 
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The components of vorticity in the r, 0, and ¢ directions are 


é= : cosec 6 cee (w sin @)| 





. r \eh 0 \? 
7 = , ma =. cosec 6 no 
r or r od 


a low 1lé 


f= 2 _- < (wn), 


respectively. In the neighbourhood of the singularity these can be expanded 


in the series af 63 

0 4 Jf 

é= me et) ts 
x) 

eee AS 
7 "adie 
63 
" 
(= (2t4g— tg) + U5 — ° 


Thus in the neighbourhood of the origin the dominant component of the 
vorticity is in the @ direction, its actual magnitude being controlled by 1. 
The flow in the neighbourhood of the singularity can therefore be regarded 
as that due to a pencil of vortex sheets of varying strength. 

A clearer understanding of the immediate neighbourhood of the singu- 
larity may be gained by use of a fixed frame of reference, and we now write 
out the expansions in terms of Cartesian coordinates. The z-axis is taken 
in the direction of the singular generator, the z-axis along the outward 
normal to this line through the vertex of the cone, and the y-axis is chosen 
to complete a right-handed system. The velocity components are u, v, w 
in the x, y, z directions respectively, and W* is the limiting value of the 
velocity attained when p = 0. 

We find that 


9 


ut = {(y—1)/6po y2wo(0)}{W*2—wR(0)}(4p,2—p' y) (x2 +y2) +. 
» = {(y—1)/6py y2wo(O)}{W*2—w3(0)}(4pyy + py 2\(v2+y) +. 
W = Wy(9)+- 309 v(y—1){ W*?— uz (0) wo g(a? +-y?)/2?-+..., 

P = Pot pal (u?+-y?)/z??-+..., 


where @ = tan-1(y/x), wo(@) = dw,/d0, (d/d6) (a2 py/wo) = 16u2 p,/wo, and 


2  19)_9 1 
w(0) = W* t -(1 _ a Jexp (ey °, 


y *2 


J 


Coy 














The 


and 


The 
sion 
P 
of t 
var’ 
the 
ang 
or § 
esti 
nel 
an 
bet 


pel 








anded 


f the 
A Un. 


ded 


ingu- 
write 
aken 
ward 
osen 
vv. Ww 


the 











A VORTICAL SINGULARITY IN CONICAL FLOW 
The dominant components of the vorticity are 


, 
Wot , 

= o 1 °°°9 

es ial 





and 1) - este epee “ee whic 


The properties of the singularity can again be deduced from these expan- 
sions. 

Probably the most significant property of the singularity for the solution 
of the yawed-cone problem is that, in its neighbourhood, all the dependent 
variables can ultimately be expressed in terms of a single unknown, namely 
the entropy distribution along the shock. This suggests that, for small 
angles of yaw at least, the problem can be reduced to one involving this 
or an equivalent unknown only. In the more general case a connexion is 
established between the shape of the conical shock and the flow in the 
neighbourhood of the singularity. This should be of value in developing 
an iterative method of integration of the equations of motion step by step 
between the singularity and the shock. 

Acknowledgement is made to the Chief Scientist, Ministry of Supply, for 
permission to publish this paper. 
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SUMMARY 

The exact ordinary differential equations of von Karman for the flow due to a 
rotating disk of infinite radius are integrated for the case of uniform suction through 
the disk. In the analysis a suction parameter a is introduced, where a,/(va) is the 
velocity of suction, v being the kinematic viscosity and w the angular velocity of the 
disk. Fora 1 the equations are integrated numerically, but for higher values of a 
a series solution in descending powers of a is obtained. 

The magnitude of the radial component of flow is found to decrease rapidly as the 
suction increases, while at the disk the derivative of the tangential component—with 
respect to distance from the disk—increases. If the ratio of distance from disk to 
displacement thickness is used as the dimensionless independent variable, the change 
of the radial component with suction is seen to occur mainly in the velocity scale, 
with little change of shape, while the tangential component of flow changes very little. 


Symbols 

a suction parameter. C 2(w/v)?. 

r, D, z cylindrical polar coordi- 7 at. 
nates. 5* displacement thickness. 

u,v,w corresponding velocity 6 momentum thickness. 
components H* 68*/6. 

F,G,H_ reduced velocity com- H,, limiting value of H(¢) at in- 
ponents. finity. 

w angular velocity. A angle of yaw in flow. 

v kinematic viscosity. Ax limiting value of A at infinity. 

p density. ph z/d*. 

p pressure. Fax Maximum value of F. 

P reduced pressure. M | 


1. Introduction 

It has been shown by von Karman (1) that the equations of steady flow 
of a viscous, incompressible fluid due to an infinite rotating disk can be 
reduced to a set of ordinary differential equations. He solved them by an 
approximate integral method, while later Cochran (2) integrated the equa- 
tions numerically for the case of zero normal velocity at the disk. The 


[Quart. Journ. Mech. and Applied Math., Vol. VII, Pt. 4 (1954)] 
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inverse problem, the case of an infinite plane at rest in fluid rotating with 
wniform angular velocity at large distances from the plane, has been solved 
numerically by Bédewadt (3) (see also Schlichting (4), p. 148). The 
slutions mentioned above are particular cases of a general class of rota- 
tionally symmetric flows which has been discussed by Batchelor (5). 

The simplified Navier-Stokes equations derived by von Karman are 
:pplicable also to the case when there is a uniform flow through the surface 
of the disk, as has been pointed out by Batchelor (5). In the present paper 
the case of suction is examined, but the case of blowing through the disk 
is omitted. For low values of the suction parameter a the solution may 
be obtained numerically, but for higher values of a a series solution in 
descending powers of the parameter is derived. In the present paper a 
numerical solution is obtained for one particular value of the suction 
parameter (a 1), while for the other values considered the series is used. 
The results are discussed in section 5. 


2. Statement of the problem 


Cylindrical polar coordinates are used, r being the radial distance from 
the axis of rotation of the disk, ¢ the polar angle in the direction of rotation, 
and z the normal distance from the disk. The respective velocity com- 
ponents are wu, v, w, while w is the angular velocity. In addition, » denotes 
pressure, p density, and v kinematic viscosity. Then the Navier-Stokes 
equations of motion (4) in these coordinates are 


- , p 
cu veu ousé 1 Op . u 2 ow 
u t + Ww = = — i +v{(V¥u—————], (2.1) 
er  réd dz fr p or rr? Od 
cv v av cv ww l op ; 20u v 
u } + ww i. — Ad + v{ V20 — —— = <=], (2.2) 
cr’ red Oz r pr ¢ r= 0 re 
ou v Ow ow le¢ ‘ 
u 4 w Pa vV2u (2.3) 
co? } od C2 p cz 
ov . Ou 
(ru) 4+—-—4+— =0, (2.4) 
r oT r oh C2 
" oc l ¢ l oe - 
where i oT —\ an) oo (2.5) 
Or? ror r= dp* a2? 
With von Kaérman’s substitutions, namely 
u rwF(C), v = rwG(C), wo (veo) H(¢) (2.6) 
, 2. 


‘menany 


p = prwP(6), C = (w/v)tz 
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equations (2.1) to (2.4) become 
F?—G@°+HF’ = F’ ) 
2FG+HG' = G" } 
oF+H’=0 | 
—HH'—2F’ = P’, (2.8 

where primes denote differentiation with respect to . 

Since H is independent of 7, these equations are applicable to the cag 
of uniform flow (w,) through the disk, H then taking a constant non-zem 


value at € = 0. The boundary conditions for uniform suction through the 
rotating disk are 


F=0, G=1 H=-—a atl=0) 
F=G@=0 at{=o f 


’ (2.7 


where a is a positive constant. 

In suction problems it is usual to introduce a suction parameter w, R/U, 
where w», is the velocity of suction, U, a representative velocity, and R the 
Reynolds number. In our case, since 


Wy = ar/(vw), Uy = ra, R= ?rw/v, (2.10) 
a is the appropriate suction parameter. 
Equation (2.8) is immediately integrable and yields 
P = constant —}H?—2F, (2.11) 
while equations (2.7) have to be integrated subject to the five conditions 
(2.9). 
For small values of a a numerical solution may be obtained, the value 


a = | being selected for illustration in section 4. A solution valid for fairly 
high values of a, of order 2 or greater, is developed in section 3. 


3. Solution in inverse powers of the suction parameter 
Let us first of all examine the form which the solution takes for very 
large values of a. On physical grounds, we expect H to be nearly constant 
for large values of a, and therefore put H = —a, with its consequences, 
H’ = F = 0, in equations (2.7). 
To a first approximation, then, these equations become 
FP’ = —aF’—@, (3.1) 
G” = —a’, (3.2) 


the boundary conditions being the relevant ones of (2.9). The solution of 
(3.2) which satisfies (2.9) is 


G = e-%, (3.3) 
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\fter substituting (3.3) in (3.1), we obtain an equation for F, the solution 
f which is 
’ ] , 
I (¢ -aG 


») 


—e-2a6), 
2a* 


(3.4) 
[he zero value of F assumed initially is seen to be an approximation to 
3.4) for large a. 

The functions (3.3) and (3.4), together with H a, may be regarded 
3 a first approximation to the solution of (2.7) for large a. The form of 
this solution suggests two things: (1) that a solution may be obtained in 
lescending powers of a; and (2) that af is more suitable than € as the 
ndependent variable. 


We therefore define ny =al 


ind transform equations (2.3) to 


a2F” — F?+aF’H—G 
GQ” = 2FG+aG'H (3.6) 
0 = 2F+aH’ J 


ere primes denote differentiation with respect to n (and not £) in the remainder 


f this s¢ ction only. 


\ solution is assumed of the form 
x 
H a+ > a*H,n) 
r=0 
F = Sa~K(n) 
G = > a'G,(n) 
\fter substituting these series in (3.6) and equating successive powers 
ia to zero, we obtain the following set of equations: 
Fo T Fo = 0, 
Fi+F; Fo Ah, 
7 , 9) 2 si sa 
F3+F, F5—Go+F, 4,+F; A: 
Cont Po, Fri» —Ga_1t+2 > (FE An-2--—G, Gen-e-r) + 


2n—1 : 
> FA, (n > 2), 
0 


2n-1-r 


r 
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n—1 2n 
FonsitFon+1 = 2 > (F. Fin-1-,—G, G,.. 1-1) + b |. is > ), 
r=0 r=0 
GQ me >= =: @, 
Gi+ _* = H, G) 
Gr +G, = Jf Gi-o-++ 3 H,G,-1-, (n> 2), 
r 
x% 
H, = —2F,-, (n> 1). (3.8) 


The relevant boundary conditions are 


F, = 0, G, = i, G, = 9, H,=0 atyn=0 ) 


(3.9) 
where y= 0, 1, 2. 3, 4...., %,... 
uw = 1, 2,3, 4,..., 2, 


After the above equations have been solved in succession, we obtain the 
following series for F, G, and H: 


1 
P= 5,3 ee) + 


+ alld Be + be 2-4 de he] + 
+ al(as7°+ 376 + 3aseu le? + (—49?— 88 — 392 e-2n + 
+(—*n-} 73 )e~39 +- - (457) — ges )e-*9-+ shh 5 — rege 4 
1 ‘ 
+O|—), (3.10) 
a 
. 1 
G = e+ —[he —17—4ne-7— 3, e-3n|+ 
at 


[(39°-+539— ifn )e-” + §(1 + n)e27+ gye-47— gh e-7] + 
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] ‘ 
B= —a4 4 (—J4e-—Je-+ 
I 201 ( 1 95 )¢ ‘ n+(4 +38)e-27+4-1 »>—3n___1 e~47| | 
7288 27 — 7 2 a4) 12° 288 a 
21023 (12-391 51829\p—7 _| ( 1,,2__119 145\,—2n 
“iil — {3980 + (89° +2887 +i7980)e "+ (— 39° —“9g' — 138) + 


i‘. 209 31 1 1 \p—4n_1. 18 9-5 7__p—67, 
r(—$7y 864 )¢ 1+ (7477 —st6)€ 1+ -57600 1—si890€ 1+ 


qs 


| o( ) (3.12) 


The numerical results for particular values of a as given by these formulae 


are discussed in section 5. 


4, Numerical solution for small values of the suction parameter 

For zero suction, Cochran (2) integrated equations (2.7) by the Adams-— 
Bashforth method (6). In order to bridge the gap between this solution 
fora = 0 and that of the previous section for a of order 2, equations (2.7) 
have been solved for a = 1 by Cochran’s method. The equations are 
first of all converted to a set of five first-order non-linear equations, which 
have to satisfy the five conditions (2.9). The solution near to ¢ = 0 can 
then be obtained in the form of series containing two arbitrary constants 
a, and b,: 
l 7 

F = a,{+-..., 

G = 1+6,f+..., (4.1) 
H = —a—a, (?+.... 


For particular values of a, and 6, the solution is then developed by forward 
integration in the Adams—Bashforth manner. 
There is also an asymptotic solution (2) for large ¢ of the form 
F ~ Ae-S-+L..., 
G~ Be-t+..., (4.2) 
H ~ —c+(2A/c)e-S-+-..., 


where A, B, and c are constants, and —c represents the limiting value of H 
for large £. The constants a,, b,, A, B, ¢ are chosen by trial and error so 
that the solution obtained by forward integration links properly with the 


asymptotic solution (4.2). 
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For the case of a = 1 the following values of the constants are applicable: 
ay 0-389, b, —1-175, A = 0-334, B= 1-034, c = 1-259, (4.3) 


The method gives the functions F, G, H, F’, and G’ directly in numerical 
form. The results for F and G are given graphically and discussed later, 

It was mentioned in the introduction that von Karman devised an 
integral (or momentum) method of solving equations (2.7). Formulae 
analogous to his can be used to check the above calculations for a = |, 
By integrating the first equation of (2.7) by parts and using the third with 
the conditions (2.9), we obtain 


x 


F’(0) = a, = [ (G2—3F*) dé. (4.4) 


0 


And, similarly, we derive 


G'(0) = b, = —a— | 4FG@ dt. (4.5) 
0 
For the case a = | integration of (G?—3F?) by Simpson’s rule yields 
a, = 03897, showing a difference of 0-18 per cent. compared with the 
value of (4.3). Integration of 4FG gives b, = —1-174, again showing a 
very small difference compared with the value of (4.3). 


5. The nature of the fluid motion 


A rotating disk acts as a centrifugal fan, the radial flow being balanced 
by an induced axial flow towards the rotating disk. When suction is applied 
the radial flow is decreased, while the axial flow at infinity towards the 
disk is larger. The effect of suction on the tangential component is much 
as for two-dimensional flow, the boundary layer being thinned. 

For a range of values of a, the functions F and G are plotted against ( 
in Figs. 1 and 2, Cochran’s curves also being given for comparison. For 
a = 1 the curves are those obtained numerically by the method of section 4. 
It may be mentioned here that (3.10) gives F accurately to four decimal 
places at a = 2, while (3.11) gives @ accurately to three decimal places. 
The considerable reduction in the three-dimensional character of the flow 
is shown in Fig. 1. The radial component shows a rapid decrease in magni- 
tude with increase of suction, the reduction compared with a = 0 being 
more than 50 per cent. for a = 1 and even larger for higher values of @. 
The maximum values of F are given below in Table 2. The thinning of the 
boundary layer is shown by both Fig. 1 and Fig. 2. 
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Standard boundary-layer lengths and parameters 

A more precise measure of the reduction in boundary-layer thickness may 
be obtained by consideration of a displacement thickness. In three-dimen. 
sional boundary-layer theory it is the practice to define two ‘displacement 
thicknesses’, corresponding to two perpendicular directions (7). In our case 
the obvious directions are along the radius and tangential to it; but since 
the radial flow is zero both at the disk and at infinity, the radial “displace. 
ment thickness’ is meaningless. The tangential component of flow gives 


5* —_— (v wr) dz, (5.1) 
in the notation of Fig. 1. With the transformation (2.6) this becomes 
8* = (v/w)' [ G dé. (5.2) 
0 


Similarly, we have a momentum thickness 


x 


8 = (v/w)! | @1—G@) dé. (5.3) 
0 
We define also H* = $*/6, (5.4) 


the superscript being introduced to avoid a clash with the notation for the 
flow function H(€). 
Within the range of validity of (3.11), (5.2) and (5.3) can be integrated 


to give ~ 
] 4 589 
8*(w/v)t = —[1— = oo 5.5) 
(w/») ‘\ 9a4* © 720a8 ) si 
ae ee | ae 5.6) 
(w, v)} = oa —_— i2a4 ose ls (5. 
. , . l i 
while the ratio yields H* — 2(1 i eee § (5.7) 
' 36a4 


For a = 0 and 1, the @ function can be integrated numerically to give 
5*, 6, and H*. The results are given in Table 1. 


TABLE 1 














a | o I 2 i = | @ x 
j re Cw Oe eee = : 
S ! | 
5* (w/v)? 1271 | o-811 | 0°488 | 0331 | 0-250 | o 
1\t 
O(w/v)? 0°599 | O-gor | 0°244 | 0166 | o-125 | oO 
2°0 2 


H* 2°122 | 2°022 
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2°000 2-000 
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Further precise properties of the fluid motion 

The axial flow H(¢) is not plotted, but its limiting value at infinite 
distances from the disk, H,,, is given in Table 2 for a series of values of a. 
It is noted that H,,+-a, the change of H between £ = 0 and f = o0, tends 
to zero as a tends to infinity. Within the range of a greater than 2, H,, is 


given by (3.12), namely 


H, a— ta-3+-20lg—-7—20zsg-U4+.. (5.8) 


Let us detine A as the angle between the direction of motion of the disk and 
the stream line projected onto the plane of the disk. Then A varies with the 
distance from the disk; it is zero at the disk, but has a finite limiting value 
\, at large [ given by 


A 


x 


- lim tan-1(F'/G). (5.9) 
(+x 

For a = | the value of X,, is derivable from (4.2) with (4.3), while a similar 

procedure suffices for Cochran’s (2) case of a = 0. And within the range 

of (3.10) and (3.11) the series are applicable. These results are given in 

Table 2, where it is seen that A,, decreases rapidly with increase of suction 


x 


and eventually tends to zero. 


TABLE 2 











° I 2 3 4 L 
poe 0181 | c:080 | 0°0295| 0°0136| 0°0078/ o 
H 0886 | 1°259 | 2°057 3019 | 4°008 | o 
6 BS¢ 0°259 | 0°057 | O'oIg | 07008 | oO 
\ q2’ | 17° 54’| 6° 47’ | 3° 10’ | 1° 48’ Jo 
In many respects » = z/5* is more suitable than ¢ for comparing velocity 


profiles. Thus in Fig. 3 the function F is plotted against yu, while in Fig. 4 
the G profiles for zero and infinite suction are plotted. The maxima of F 
are at about the same value of y, and a large part of the change of F with 
suction is seen to be one of scale. To show this the F profile for large a, 
namely 
F 33° He). (5.10) 
is scaled up to lie near the F curve for a = 0; the closeness is remarkable. 
Similarly, the F curves for a = 1, 2, 3 do not differ much from suitably 
scaled versions of (5.10). In deriving (5.10) one notes from (5.5) that 7 
tends to u for large a. 

The curve for G against ~ at a = 00, which is G = exp(—y), is not very 
different from that for a = 0, only one suction curve being plotted as they 
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all lie close to that for a = oo. But the qualitative change from a = 0 to 


1 = © is similar to that for two-dimensional flow. For a Blasius distribu- 


tion the parameter H* changes from about 2-6 for zero suction to 2 for the 


asymptotic suction profile, whereas here H* changes from 2-122 to 2. This 
small change of * is seen to be in agreement with the closeness of the 
curves of Fig. 4, since the value of H*—2 is a measure of the closeness of 
, velocity profile to the corresponding asymptotic suction profile. It is 


shown by Figs. 3 and 4 that, in the variable py, the velocity distribution 
with or without suction can be represented approximately by 


F = M(e-+—e-2¥), 


(5.11) 
G e--, 
where M _ (5.12) 


These formulae become increasingly accurate as a increases. But, of course, 
these approximate formulae do not hold as closely for derivatives of F and 
Gas they do for the actual functions F and G. 
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SUMMARY 

ey P ; ; / 

The non-linear hydrodynamical equations for the steady axisymmetric flow of 
perfect fluid have been reduced to one linear differential equation for Stokes’s 
stream function. It is shown that uniform motion of bodies of revolution along 
the axis of the rotating liquid can be discussed on the basis of this equation. 


1, Introduction 
THE motion of solids in rotating liquid was studied in a series of papers 
published in the decade 1916-26 by G. I. Taylor, J. Proudman, and §. F. 
Grace. With the exception of Taylor, these writers have usually linearized 
the equations of motion by assuming that the liquid rotates with a uniform 
angular velocity about the z-axis, which is also the axis of symmetry of the 
body of revolution, and that the body moves with a small uniform velocity 
along the z-axis. The linearized equations form the basis of a series of 
interesting investigations recently carried out by Gortler (1), Morgan (2), 
and Stewartson (3). In one case, however, Taylor (4) obtained a family 
of solutions of the non-linear equations when a sphere moves with a uniform 
velocity U along the axis of the rotating liquid. He reduces the flow toa 
steady motion by superposing a velocity —U on the whole system and 
obtains an infinity of solutions satisfying the kinematical boundary con- 
ditions, all of which have the property that the relative velocity vanishes 
at infinity. 

In this paper we show that the non-linear equations for steady motion 
of a body of revolution along the axis of a rotating liquid can be reduced to 
a single linear equation, which may be solved for any body of revolution. 


2. Equation of motion 
Consider a body of revolution moving with a uniform velocity U along 
its axis of symmetry and the liquid rotating with a uniform angular velocity 
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Q about this axis. The flow can be reduced to a steady motion by super- 
posing a velocity U on the whole system. Hence the body comes to rest 
and the liquid moves with a velocity —U at infinity along the axis of 
symmetry of the body and also rotates with a uniform angular velocity Q 
about the axis. Taking the origin at the centre of the body and the z-axis 
as the axis of symmetry, the equations of steady motion are 


grad(! +f) —Vxcurl V 0, (1) 
2 p 
div V = 0, (2) 


where p is the pressure, p is the density of the liquid and V the velocity. 
Since the motion is axisymmetric, we follow Goldstein (5) in introducing 
Stokes’s stream function, % and replace (1), (2) by their equivalents in 
orthogonal curvilinear coordinates a, B, y. We find 


l r us l Cus 





€u : v= —— : haw = y, 
hahz op hyhs ex : . 
(3) 
; l oéx l ox P 
é A een ; hal = —D*f, 
; hahs op 7 hy hs tale " id 
oan D2 h, [é ( h, @° re h, h 
hy hg|@a\hgh, ex OB\h, hs ep 
The third (y) component of the equation for the vorticity is 
2y O(y,h3) Ob, D*b) | 2D A(x, hg) _9 (4) 
h. O(a, B) (a, B) h O(a, B) 7 
The third component of the equation of motion is 
Ol, x) __ 9, (5) 


O(a, B) 


In these equations (a,8) are the general orthogonal coordinates in the 
meridian plane and y is the azimuthal angle; h, is the distance from the 
axis of revolution; (w, v, w) are the components of the velocity vector V 
in the directions of a, 8, and y increasing: (€, , £) are the components of 
the vorticity vector w. 





From (5) we have x = f(p). 


4 


In view of the boundary condition and the condition at infinity on w we 
take 


f(b) = —hf, 
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where k is a constant. Then, 


w = kib/hs; (6) 
) p2y fe 
and (4) reduces to ali ad = (7) 


where H is an arbitrary function. If the liquid rotates with a uniform 
angular velocity Q at infinity, we have 


[w).. = hg Q, 
i.e. [eb]. = —AZQ/k, (9) 


where the symbol [y],, stands for the value of % at infinity. 
Using (9), equation (8) becomes (cf. 6) 


(D?+k®)b = —kOR, (10) 
whence = — Oh? k-4 Ay, + Bobo, (11) 
k= 20/0, 


where ys, and #, are the two solutions of 
(D?+k?)b = 0 (12) 


since D*h3 = 0, as is shown in the Appendix. 
In any specific problem, h, being known, it is only necessary to solve 
(12) in a suitably chosen system of orthogonal coordinates. 


Case of a sphere 
When the body of revolution is a sphere, we use spherical polar coordinates 
(r,6,¢). Putting «a =r, B = 0,h, = 1,h, = 1, hg = rsin8@, (12) becomes 


eu sind é 
er r2 dé 


(cosee 0m) ao ky = 0. (13) 
00 
Then the solution (11) is of the form 
Y = —}U r*sin?0-+ (kr)! S sin 6 P},(cos 6)[A,, J,,,,(kr)+B,, J_,,,(kr)]. (14) 
n=1 * 


This agrees with the result of Long, who obtained and solved (10) in 
cylindrical and spherical coordinates. Taylor’s (4) solution can be obtained 
from (14) by taking the first term of the series corresponding to n = 1. 
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HYDRODYNAMICAL EQUATIONS 


Case of a spheroid and a circular disk 

Equation (12) has been solved by Fadnis (7) in spheroidal coordinates 

terms of spheroidal wave functions. He has discussed the motion of 
oblate and prolate spheroids and deduced the motion of a circular disk as 
a limiting case. 

The authors wish to thank Sir Geoffrey Taylor for informing them of 
the work of Dr. Long and for clearing up some doubtful points. 


APPENDIX 


Consider the transformation 
r+iz = f(a+ iB), 
whet nd z are the cylindrical coordinates. 
If .+i8, and @ be the conjugate of w, we find 
j I f’(w)+f'(@)], hy Af (w) -f'(w)), hs 3 f(w) +f (@)] 
t(w) f(a {l < C , ae 
pg =} LENO (2 4 2)ey 0) y10y- 
2 t**(w) t“(a) Cw Ccwi~ j 
c c A ee 
oni ——|[f’(w)- f (®))| 0. 
Cw Ow 
re D?) 0 


,EFERENCES 
H. GOrtTLER, Z. angew. Math. Mech. 24, Nos. 5 and 6 (1944). 
G. W. Moracan, Proc. Roy. Soc. A, 206 (1951), 108. 
TEWARTSON, Quart. J. Mech. and Applied Math. 4 (1953), 141. 
G.I. Taytor, ‘The motion of a sphere in a rotating liquid’, Proc. Roy. Soc. A, 102 


- wh 
} 
I 


1922). 180. 


STEIN, Modern Developments in Fluid Dynamics (Oxford, 1950), vol. 1, 


on 
I 


6. Ros R. Lone, Technical Report No. 1, Johns Hopkins University (Dec. 1952). 


ADNIS, unpublished. 














ON THE PROBLEM OF DIFFUSION FROM AN 
INSTANTANEOUS POINT SOURCE RELEASED 
AT GROUND LEVEL INTO A TURBULENT 
ATMOSPHERE 
By D. R. DAVIES (The University, Sheffield) 
[Received 14 May 1953; revised 25 June 1953] 


SUMMARY 


In this paper it is shown that O. G. Sutton’s prediction formula, describing 
diffusion from an instantaneous point source released at ground level into a turbulent 
atmosphere, is a particular solution of the three-dimensional equation of mean 
atmospheric diffusion, satisfying the appropriate continuity and boundary con- 
ditions. The eddy diffusivity coefficients employed in the diffusion equation are 
expressed as power law forms of the time of diffusion (measured from the instant 
of generation of the cloud), the indices and proportionality constants in the power 
laws being specified by an application of Sutton’s extension of a theorem on diffusion 
due to G. I. Taylor. 


1. Introduction 

FoRMULAE describing average diffusion from instantaneous point sources, 
released at ground level, have been obtained by O. G. Sutton (1, Chapter 8): 
these satisfy the equation of continuity imposed, the boundary conditions, 
and also the condition prescribed by Sutton’s extension (1) of a diffusion 
theorem due to G. I. Taylor (2). These formulae have been shown (1, 
Chapter 8) to be satisfactory for practical purposes of prediction in con- 
ditions of neutral atmospheric stability. Sutton, however, has indicated 
(1) that the diffusion expressions are not unique, and they have r-t been 
shown to be solutions of an appropriate diffusion equation. 

In the present paper the problem is approached of obtaining a solution 
of the three-dimensional equation of atmospheric diffusion subject to the 
conditions describing the initiation of an instantaneous point source at 
ground level, the parameters involved in the forms assumed for the eddy 
diffusivities being determined so that Sutton’s extension of G. I. Taylor's 
diffusion theorem is satisfied. 


2. The three-dimensional equation of atmospheric diffusion and 
boundary conditions describing the initiation of an instantaneous 
point source at ground level 
Referred to fixed rectangular axes, Ox at ground level in the mean down- 

wind direction, Oy in the lateral direction, and Oz vertically upwards, the 


(Quart. Journ. Mech. and Applied Math., Vol. VII, Pt. 4 (1954)] 
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DIFFUSION FROM AN INSTANTANEOUS POINT SOURCE 


} (1) 


where y is the average concentration of air-borne particles, U the mean 
wind speed at a level z above ground, and K,, K,, K, are the diffusivities 
in the x-, y-, 2-directions. This equation is taken to describe the stochastic 


full mean diffusion equation is 


Ox , yy OX c (Ke) + (K+ (ie: 
Cx 


at Ox “Oa) © dy\ © dy) © "| 


N ioe 


Oz 


average diffusion over a large number of separate experiments involving 
the release in each experiment of a single instantaneous point source, in 
identical conditions of neutral atmospheric stability (apart, of course, from 
the random fluctuations inherent in the turbulence problem). The dif- 
fusivity coefficients K,, K,, K, are regarded as statistical quantities ex- 
pressed in terms of the relevant variables and mean parameters of the flow. 

The boundary conditions necessary to describe the release of an instan- 


taneous point source at the origin are: 


x>0 ast-0, atr=y=2z=0; (2) 
x20 ast>0, atx>0,y>0,z>0; (3) 
x > 0 as > 00, Y¥> 0, 2> 00; (4) 


and the condition of continuity is that the total mass of air-borne particles 
should be constant and equal to Y, the mass emitted instantaneously at 
the origin at ¢ = 0. 


3. Solution of the diffusion equation in terms of power law 

diffusivity forms 

The physical notion that progressively larger eddies dominate the mean 
outward diffusion of the cloud as the time of diffusion of the cloud increases 
suggests that the coefficients of eddy diffusivity K,, K,, K, are to be regarded 
as functions of time, measured from the instant of generation of the cloud. 
We now assume that these coefficients are expressible as power-law forms 
of the time, and we write 


K.=<4,®, (5) 
K, a y t . ? (6) 
K, = a,#*, (7) 


where the index « and the proportionality factors a,, a,, a, are constant for 
a given turbulent flow field and their values are to be determined. 

It has been shown by Sutton (1) that the errors involved in prediction 
formulae are small, if we assume that the mean wind speed U is independent 
of height z in considering the term U @x/éxr on the left-hand side of the 
diffusion equation (1): the observed variation of U over most of the normal 
height of cloud in practical cases of importance is fairly small when neutral 
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stability conditions are considered. In this paper we consequently assume 
U = U, (constant) in equation (1), where U, is the observed value of the 
mean wind speed at a convenient reference height z,, which is taken to be 
2m in practice. 

Substitution of equations (5), (6), and (7) into (1) leads to the equation 


| < 


C » © G o* co” 
X41 ox = t/a, - Xia, —A+4 _ x : (8) 

Ca Ca? ” oy? 62? 
We now introduce the new variable 


. 


€é=2-—U,t. (9) 
Extending the method of a previous paper (3), we write 
Evt-A = U, yt A — y, 2t-A w, (10) 
and seek solutions of the form 
x = At™f(u, v, w), (11) 
where A and n, are important diffusion parameters, whose values are to be 
determined. The combination éy/ét+ U, éy/éx, which appears on the left- 
xX nad 
hand side of equation (8), is found to transform to the expression 
C C of 
Afn, t”°-1— AXt™o-! u F . v ao w I . 
ou ov ow 
and it is interesting to note that this is valid even if U, is a function of z. 
The resulting differential equation is not separable, however, unless we 
1 ; 
assume U, to be a constant, owing to difficulties in dealing with the term 
involving @y/éz. 
Substitution of (10) and (11) into equation (8) leads to the equation 


a | (t-te) of) oii. 
u 


C ul 


‘ C C C , 
fngt" 1__)gro-1l y : f Ly Sf of A, P't*+No-2Ap)y (1 1/p) 
| OU ov cow ; 


! a, GX +o 2A/q)py(1 1/q) vd fa 1/q) ef | M a. r2t(*+No 2A/r)94,(1 l/r) c uf AN 
; ov ov 7 Cw cw 


(12) 
We now choose values of the parameters so that the ¢ variable is cancelled 
out in equation (12) and we obtain 


9 9 9 
2A _ 2A ; 2A l Aw, (13) 
Pp q r 
We note that the continuity equation is 
4 ( ( ( x dédydz = Q, for all values of f¢, (14) 
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ind this gives 


4 Atrot3Ap))-3 ( [ f (u,v, w)UP-) dudvdw = Q; 


0 0 0 


hence No 3A/p —3(1+4)/2, (15) 


using equation (13). 


The differential equation to be solved is now 


| , at 
Ny J—A\U 2 4@ J w A a, p?ut-lp) be | yt-a oJ | + 
: cu ov COW , cu eu} 
a, p2va—-lp) © } ya ip) FV, a, p2w%t-Up) ef yan BF), (16) 
v | ev} : ow | éw| 


Writing f = F(u) G(v) H(w), we find equation (16) is separable, and we 
btain the following normal solutions, in exponential forms, 








' (1+-a)u2!p 1 +- a) GeO +) » 
Pw) = exp —O-tSIe _ geg/ te ay 
| 40, } | 40, J 
(] 2/p 1+ 2¢—(L+a) 
(a) exp! ai | exp| —E TY es (18) 
| $a, 20, J 
(id x) wr! | (1-4 x)2%e—“+e) 
H(w) = exp — exp) - —-——_ |. (19 
I fa. | P) 4a, 
\ solution of equation (8) is then 
- \ £3 ye 22 
y = At-M+artexp)_ ETS ZY (20) 
| 4@+%]@, a, a,j} 
where A is found from equation (14) to be given by 
A (l+a)iQ 4ri(a,a,,a.)'. (21) 


The parameters a, @,, a@,, a, are now determined by applying O. G. 
Sutton’s extension of G. I. Taylor’s diffusion theorem. We evaluate ex- 
pressions for the standard deviations (o¢, o,, o,), from the cloud centre, of 
the particles lying, at time ¢, along the three lines which are parallel to 
Ux, Oy, Oz respectively and pass through the ground-level centre of the 
loud. The result in the mean wind direction Oz is given by 


a = | x2 ae/ | x dé, (22) 
0 0 


with y and z taken to be zero in the y expression given by equations (20) and 
21). Hence 


: 2 
of = 24_(1 17+% noting that K, oi. (23) 
é a, X : ; 9 dt é 
3ut o2 is given also by Sutton’s result (1) that 
o 102 U0 m)¢g2i(l+m) | (24) 


92.28 Hh 
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where C; is Sutton’s generalized coefficient of diffusion (1) in the down. 


wind direction and m isa turbulence parameter calculated from the observed 
mean wind velocity profile. Equating results (23) and (24), we obtain 
x = (1—m)/(1+m), (25) 

' - 2 [j2(1+m) 
and a, = }(1+ x) OF L ere, (26) 

By considering now the standard deviations about the centre in the Oy 
and Oz directions, we find that expressions similar to that found for a, are 
obtained for a, and a, with C, and C, respectively replacing C;. 

With the values found in this way for «, a,, a,, and a,, the final form of 
the instantaneous point-source solution is then seen to be 


20 nal 1 E / . al ‘i 


= — ~ exp) — —______ | —.+ ~—.+-]}, 2] 
mC: C, C(U, t)a+™ P| (l £4] CF" ( : "@ | 





which satisfies the normal form 


ol i i 
' ed aleatoat cal] 


y % 
20 (87°)'o¢ 0,0, 





This is precisely Sutton’s prediction formula (1) for an instantaneous 
point source, and was derived originally by Sutton to satisfy the equation 
of continuity imposed, the boundary conditions, and the condition pre- 
scribed by his extension of G. I. Taylor’s diffusion theorem. We have now 
shown that this formula also satisfies the three-dimensional equation of 
atmospheric diffusion, with U, taken to be a constant and the diffusivities 
in the x-, y-, and z-directions dependent on a power-law form of the time 
from the instant of generation of the cloud. As far as the author is aware, 
no reliable experimental results are available to test the diffusion formula 
(27) directly, but, as shown by Sutton (1), a solution describing diffusion 
from a continuous point source can be obtained by integrating equation 
(27) with respect to time from zero to infinity. The order of agreement of 
the ensuing theoretical results with the results of diffusion experiments in 
neutral stability conditions is found (1) to be good; e.g. the predicted 
dimensions of cloud width and height at 100m. from the source are 40m. 
and 13m. compared with the observed values of 35 m. and 10-5 m. respec- 
tively, these being time means over periods of at least 3 minutes’ duration. 


4. Conclusions 

In this paper it is shown that O. G. Sutton’s prediction formula, describing 
diffusion from an instantaneous point source released at ground level into 
a turbulent atmosphere, may be obtained as a particular solution of the 
basic three-dimensional equation of atmospheric diffusion (with U, taken 
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to be a constant on the left-hand side). The eddy diffusivity coefficients 
employed in the equation are assumed to depend on the time measured 
from the instant of generation of the cloud, and are written in the form 


(1+? fl—m)( m). K 1(] | an) CZ UR + gm +m) 


9 
1/] yCzT y 4 1 


gv 
K. 1(] | x) C2 U2/A mf m)(1+m) | 

the proportionality constants and indices having been determined by an 

application of Sutton’s extension of a theorem on diffusion due to 


G. I. Taylor. 
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THE PROBLEM OF DIFFUSION INTO’ A TURBULENT 
BOUNDARY LAYER FROM A PLANE AREA SOURCE, 
BOUNDED BY TWO STRAIGHT PERPENDICULAR 
EDGES 
By D. R. DAVIES (The University, Sheffield) 
[Received 7 October 1953; revised 11 December 1953] 


SUMMARY 
In this note a solution of the equation of turbulent diffusion in the atmospher 
is obtained, subject to the boundary conditions prescribed by a plane area source 
in the shape of an infinite quadrant, when the surrounding surface is assumed to be. 
impervious to the transfer of particles and the mean wind speed is directed parallel 
to one of the edges of the quadrant. 


In this note we consider the problem of obtaining prediction formulae for 
the distribution of the rate of evaporation and of the distribution of mean 
concentration of vapour associated with a plane area source in the shape 
of an infinite quadrant, on the supposition that the surrounding surface is 
impervious to the transfer of particles and that the mean wind speed is 
directed parallel to one of the edges of the quadrant. We assume that the 
vertical diffusivity A, and the lateral diffusivity K, are respectively pro- 
portional to z1~” and to 2”, where z denotes height above the surface and 
m is the parameter involved in the mean wind-speed power law of the form 
u/u, = (z/z,)", u, being the mean wind speed at height z,. It has been 
shown empirically by the author (1) that a diffusion model expressed in 
terms of these power-law forms may be used to describe adequately the 
distribution of mean concentration over an evaporating area of parabolic 
shape. This model has also been shown by the author to lead to theoretical 
net rates of evaporation in good agreement with experiment for small, 
saturated, plane, parabolic-shaped areas in a wind-tunnel boundary layer. 

If U denotes the mean wind speed at height z above the surface in the 
positive sense of the z-axis of rectangular axes O(x,y,z), and K,, K,, A; 
denote the coefficients of diffusivity in the x, y, z directions respectively, 
the equation satisfied by the mean concentration y(x,y,2) of vapour in 
problems involving uniform, continuous diffusion is 


ox ee: i O/,- 0x 1) 
‘ge. K, x) +i 14, x). (1) 
Cx Cx Cx cy cy C2\ C2 
We neglect K,, and using the notation of previous work (1) we write 
U U, = 2 Zz)", k, a, U} Nyl-m and K, a, U} nym 


[Quart. Journ. Mech. and Applied Math., Vol. VII, Pt. 4 (1954)] 
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where ” 2m/(1-+-m), and the value of a, is chosen empirically in the 
manner previously described (1). Equation (1) then becomes 
cx (“1 =#)- = tar (2 "5P (2) 
Ca \ Ut Jey oY, UF jaz Cz 
We now suppose that the infinite quadrant formed by the positive x-axis 
nd the negative y-axis is a saturated evaporating area, and that the 
remainder of the plane is impervious to the transfer of vapour. If the 
oncentration of vapour on the evaporating area is denoted by x, (the satur- 


tion value), the boundary conditions then give 


¥/xX9 > 1 x > 0. y <0, e—> 0. (3) 

K.éy/éz > 0 x > 0. y > 9, 2-0, (4) 

. 0. ~~ > ©. 2 0. (5) 

nd x > 0 2>0, 2>0. (6) 


We now write 


Ur \} » (ur\s 
r, n=(—t)y, ¢ (—2)° 2m (7) 
A, 24 (2m--1)\a, 27" 
ind introduce polar coordinates such that 
p cos 8, C psin @. (8) 
Equation (2) becomes 
Cr C"+ | 8) CX 8s CX | c*y 
“1+ -— cot ee - (9) 
cé Cp* p Cp p* C p~ c= 
where s 1/(1+2m). When s 0, the problem reduces to the heat con- 


luction case (2, p. 166). In the atmosphere s is not zero, but the methods 
leveloped for the heat-conduction problems may be extended so as to deal 
with equation (9). This equation is separable, and suitable normal solutions 
exist in the form of exponential, Bessel, and Associated Legendre functions, 
viz. exp(—a?é), p—#8J vp), and (sin @)!4-9 P21—5)(cos @); n is taken to be 
1 positive integer and a is a separation parameter. 


\ solution of equation (9) is thus 
¥ | e-*%p- J, .,(ap)(sin 6)44-9 P24—8)(cos 0) f,(a) da, (10) 


where the f, («) functions are to be determined by the boundary conditions. 


T 
y 


polar coordinates these are 
y/Xo > 1 as 0—>7, (11) 


(sin @)§ey/ce@ > 0 asd 0, (12) 
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while (sin @)°éy/@6 > a non-zero constant, multiplied by a function of » 
and €, as 6 > =z, and 


x > 0, &> 0, C>0. (13) 


We now use the formulae for P'~*)(cos @) given by Hobson (3, equations 
55 and 58), where $(1—s) and m may be non-integral and cos @ is real, and 
find that 


(sin @)§4-*) P#1—5)(eos 6) > 240-9/T(4+4+4s) as 0>0, (14) 
and (sin 0)84—-9 P11—8)(cos 0) +0 as 0—> 7. 


Moreover, as 6 > 0, 


(sin 6)§ . 


sg L(sin 0)" ®) P¥1—8)(cos 8) | > 0, 


and as 6 > 7, 


(sin @)8 . 


7p Lisin 8)!» Pit (cos 4) 


> F{—n,n+1; $+-48; 1}/{26+4P (44 3s)} 


and hence the P'~*)(cos @) functions are suitable in this problem and equa- 
tion (10) fits the boundary conditions (11) and (12); F is a hypergeometric 
function. It remains to determine /,(«) from the leading edge condition 
(13), i 


3 (sin @)§1-* P1—8)(cos 8) { p 18 J (ap) f(x) da = 1. (15) 


0 
n 6 


We first write f,,(«) = a, «~**-!, where a, is a constant, and using Weber's 
infinite integral (4, p. 391) we find that 


x 


a, 1(4n+4—}s) 





js—1 Seb A - sav 5 
a, | J,.4(xp)(xp) 8-15 da EPG inte) =b,, say, (16) 
and we are left with the equation 
l > 5,(sin @)'4-9 P#2—*)(cos 8), (17) 
n=0 


to determine the b,, coefficients. 

The orthogonal properties of the Associated Legendre functions PX? ~*(2), 
when 3(1—s) is fractional, may be easily obtained, using Hobson’s formulae 
(3). Writing cos @ = z, we find, provided }(1—s) < 1, that 


b } (1—z2?)-#a 


6) PX 8)/ ~ 2ate| f [ Pia $)(z) |? dz. 
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The final solution of our problem may then be written 





x 
y/Xo | >: (sin @)*2—-* P#1—8)(cos A)a,, | € BT, .4(u)u is-1 dy, (18) 
n 0 4 
2 ¢-—t « 1 
where 6 pé—*, anc 
1 
21+isP(Ln | : | 1s) | (1—z? i(1 s) pia 5)(z) dz 
“ 
Ai ; 7 a a aaa 
(dn | 4 —1s) [ Pi *)(z)]? dz 
1 


We note that y is a function of 8 and @ only, and that 


| e—u* By 4 (uu is—1 dy 
0 
sconvergent provided s < 1. The lateral spread of vapour at ground level, 
0, is easily obtained from (14) and (18) to be 


x x 
)3(1—8) 


tite Bae bs) 2 a, |e PS, .s(uju-'-! du, 


(3 ‘ 
™ n=0 0 
ind the distribution of rate of evaporation E in a corner of a quadrantal 


wea is give n by 


- : (7? 1 
Ei(p, £) a.l "( l (m- 1 )8p8 ly 
2 a 5 


aw it 
| 
: x 
= F{—n,n+1;4+-48; lia, [ si 
ed a | CUP (eds de. (19) 
Lat $01 sT(2 | \s) n 
n=0 sa a : ° 
0 


Experimental results are not available to verify these formulae, but it 
has been found by computing a special case (s = 0-80, m = 13) that the 
speed of convergence of the series involved and their convenience for 
computation is satisfactory. Finally we note that these formulae will give 
an approximate description of the distribution of rate of evaporation and 
of mean concentration in the neighbourhood of a leading edge corner of a 


inite, rectangular, evaporating area. 
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THEORY OF THE SQUIRREL-CAGE INDUCTION 
MACHINE DERIVED DIRECTLY FROM 
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By E. MISHKIN (Israel Institute of Technology, Haifa) 
[Received 16 July 1953] 


SUMMARY 


The current-vector locus found experimentally deviates substantially, in form, 
from the classical circle because the ordinary theory, wrongly, assumes a constant 
leakage inductance. The squirrel-cage’s preponderant skin effect causes the leakage 
inductance to vary with the slip. 

Due to the alternating teeth and slots, the teeth zones of the stator and rotor 
exhibit different mean permeabilities in the radial and peripheral directions, A 
machine model is chosen which substitutes the teeth zones by an anisotropic 
magnetic mass with equivalent permeabilities. This model, while sacrificing the 
finer details of construction, yields, through direct integration of Maxwell’s field 
equations, a unified analytical solution of the machine’s overall performance at 
variable slip. The stator current and torque are calculated as functions of the slip 
and compared with those obtained by the classical method. 

Formulae derived for the stator current at no-load and infinite slip prove to be 
very simple and suitable for the designer. 


List of symbols 


permeability of free space (47 10-* henrys/em.). 
relative permeability of teeth area. 

relative x-permeability of the anisotropic mass. 
relative y-permeability of the anisotropic mass. 
conductors’ actual conductivity (mhos/em.). 
conductivity of anisotropic mass (mhos/em.). 
linear current density (amp/cm.). 

current density (amp/cm.’). 

pole pitch (cm.). 

number of slots per pole and phase. 

number of conductors in each slot. 

slip. 

angular frequency. 

electric field intensity (volts/cm.). 

electric field intensity in the overhang zone (volts/cm.). 
magnetizing force (amp/cm.). 

magnetic flux density (webers/cm.’). 

length of stator’s mean turn (cm.). 
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core length (cm.). 


1. stator overhang conductor’s length (cm.). 


] stator current (amp). 
[, no-load current (amp). 


[,, stator current at infinite slip (amp). 


1. Introduction 

THE squirrel-cage induction machine is the most widely used industrial 
motor. Graphical methods according to the well-known circle diagram are 
often used in order to forecast the machine’s performance. In order to draw 
the circle diagram two quantities must be known: (1) the magnetizing 
current, (2) stator and rotor leakage reactances. 

The classical theory assumes a constant leakage inductance, but this 
assumption is basically wrong because of the iron saturation and the 
squirrel-cage’s preponderating skin effect. Thus these two features, which 
have not been taken into account, cause the current diagram, found experi- 
mentally, to deviate substantially fiom the classical circle. 

The aim of this paper is to find a unified analytical solution of the induc- 
tion machine’s performance at variable slip through direct integration of 
Maxwell’s field equations. By this method the whole conception of induc- 
tance may be dispensed with and full account is taken of the skin effect. 

The following theory does not take into account the hysteresis and satura- 
tion effects as saturation leads to non-linear equations and the method 
adopted uses linear systems exclusively. The results obtained could never- 
theless be corrected for these effects. 

The final formulae express directly the current as a function of the slip. 
For numerical calculations, tables of hyperbolic functions with complex 
variables often have to be used. The limit cases, ‘no load’ and ‘infinite slip’, 


yield very simple expressions (see (49), (50)). 


2. The machine model with an anistropic magnetic mass 

Fig. la shows part of a developed machine. An orthogonal system of 
coordinates is used whereby the a-axis coincides with the stator’s inner 
surface, the y-axis points towards the rotor, and the z-axis, parallel to the 
machine’s shaft, forms with the two previous ones a right cartesian system. 

We shall assume the iron in the stator’s and rotor’s yokes to be perfectly 
permeable, while, due to high induction values in the teeth zones, the per- 
meability of the iron in these zones should be assumed to be finite. We shall 
denote its relative permeability by p. 

The whole teeth area reveals a different magnetic reluctance in both x- 
and y-directions. For the purpose of the following analysis the original 
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alternating teeth and slots may be replaced by a homogeneous anigo- 
tropic magnetic mass with different relative permeabilities ,, in the x- and 
}, in the y-direction respectively (see Fig. 15). 
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Fic. la. Part of developed machine. 
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Fic. 1b. Motor with anisotropic magnetic mass. 


Actually, ,, the permeability in the x-direction, is to be represented by 
a periodic, rectangular-shaped function of x. In what follows only the 
constant term of its Fourier series is being considered; the other harmonics, 
while revealing some finer details, lead to Mathieu equations. The general 
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theory, which is being derived here, is not influenced by the higher har- 
monics. Let 
d = air gap length (cm.), 
d' = lengthened air gap (cm.) due to indented bounding iron surfaces 
(1), (2), 
d, = slot depth (cm.), 
b. = slot width (cm.), 
teeth width (cm.). 
The anisotropic mass will be the magnetic equivalent of the original teeth 
zone if 
’ 
ee oy [(0, +); ay (6, uby [{o,-+b,+ =, 
/ be d, J 
(1) 
The above relative permeabilities were determined by Ollendorff (3) by 
postulating that in the anisotropic magnetic mass a magneto-motive force 
M acting in the x- or y-direction should produce in this direction the same 
flux ® as it would have produced in the original teeth and slots zone. 
From the electrical point of view, the above hypothetical magnetic mass 
must be conductive in the z-direction only, so that the model will represent 
faithfully (although in the form of a continuous distribution) the axially 
flowing electric currents. 
It is clear that the attributed conductivity o (mhos/cm.) should be 
o bs 6, (2) 
b.+b, 
where ¢ is the machine conductors’ actual conductivity and y is the slot’s 
space factor. 
In what follows the indices 1 and 2 will be used for stator and rotor 


quantities respectively. 


3, The differential equations of the electro-magnetic field 

In what follows we are assuming a three-phase machine with a symmetri- 
cally wound stator. At present we are neglecting the end effects (stator’s 
overhang) which will be dealt with in section 6. Suppose all the stator 
conductors to be evenly distributed on the stator’s inner surface. A current- 
carrying surface with all currents flowing in the z-direction would result, 
its linear current density (A amp/em.) being composed, as is well known, 


of harmonics 1, 5, 7,..., the first harmonic being 
3v2w!l sin(7/6) ; ; 
A A — rese wl—72z 7 _ Pe \ re eiot—rr 7), (3) 


| 7sin(7 6q) } 
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where re denotes the real part of the complex quantity, j = ,/(—1), 7 is the 
pole pitch, w is the number of conductors in each slot, w is the angular 
frequency, and q the number of slots per pole and phase. 

In all that follows, re and the factor ei will be omitted, being self. 
understood. 

Actually, the conductors lie, not on the stator’s surface, but embedded 
in fairly deep slots. In lieu of the linear current density A, the current 
density vector 7 (amp/cm.?) must be introduced and, due to the winding 
method, it should be considered to be independent of y in the stator’s teeth 
region. Hence we have 


. A we 
t= — = tax OC? ™".. (4 
dy, max ) 
Maxwell’s second equation 
. 0B ; 
curl E == —jwB (5) 
Cc 


expressed explicitly in the form (E has only one component, in the z-direc- 
tion, while B has no component in this direction) 





i’. & &! 
o C C , i ; ‘ . 
— = i * jopoli, wh, H,+i, w, H,+i, 9) (6 
| 0 O E 
leads to the equations 
==) OF ] On - 
H, aap hag, (7) 


Jpg By Ox 
as B, = pou,H, and B, = pop, H,. where po is the permeability of the 
free space. 


Jpg by CY 


Maxwell’s first equation curlH = i (8) 
expressed explicitly, in the stator area, yields 
. tee 
| € c c — 
_ — |= i; 
CX CY ‘i 
|H, H, 9 
eH, eH , 
hence ¥__=—4 (9) 
Ox cy 


as the current flows in the z-direction only. 
Combining expressions (7) and (9), one obtains the differential equation 
of E in the stator’s teeth region, namely 
&E eH 


Bayz 72 Ty dy? a JLo Mr, by, i, o> ae dy. (10) 
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the air-gap the relative permeabilities u.. and yw, are both equal to one, 
n I | He I 7] 


and the current density 7 zero; hence 


“" ae ste d>y>0. (11) 
ct” oy 


In order to find the differential equation of Z in the rotor teeth zone, a new 
rotating system of axes O'(x', y’,z’) should be used, parallel to the preceding 
system but at rest with regard to the rotor. The new system is defined by 


the relations 


x’ x -<9 8)wt; y’ Y; e = £. (12) 
In the new system of coordinates the angular frequency is no longer w but 
sw where s is the slip. The required differential equation follows directly 
from (10) since Maxwell’s laws remain in force in this system but with the 


changed frequency 


Perse ths J J8@[g by Mygt23 d+d,>y>d. (13) 


4. The boundary conditions 


Fig. 2 shows the developed rotor, as seen from the stator, together with 


its end rings. The squirrel-cage bars are not shown as they are replaced 
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Fic. 2. The developed rotor (seen from the stator). 


+] . ° ° ° ° P 
by the anisotropic inductive magnetic mass (see section 2). Let r be the 
resistance of one ring per cm. length (Q/em.) and J the ring current per 
em. in the y-direction (amp/em.). The rotor being short-circuited, 
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neglecting the field at the rotor’s ends, the line integral along the contour 
shown in Fig. 2 yields 


“Cigte) ig(%)]+2 | L’rd, dx’ = ILE'(x,)— E'(x3)], (14) 


. 


Fo 


ry 


and as x, tends to x; this takes the form 


L 1; ok’ 
—— — + 2)'rd, = —I—. (15) 
Oo Cx Cx 
rT . . . . . 
The current continuity at the rotor’s ends implies that 
ol’ “7 a 
—> —— Sie (16) 
Cx 
Differentiating (15) with respect to x’, one obtains 
Ll a0, 7 ek’ 
nw et Sion, =: —f—.. (17) 
o da’? Ox"? 


The solutions of the differential equations (10), (11), (13) of the previous 
section must satisfy the following boundary conditions: 


(i) Due to the perfect permeability of the stator and rotor yokes, the 


tangential magnetic field H, must vanish along the planes y = —d, and 
y = d+d,. Hence, considering equation (7), we obtain 
(Fr) os (: 7) ont: (18) 
cy y=—dse oy y=d+d, 


(ii) Along the boundary planes y = 0 and y = d the tangential intensity 
H,, and the normal induction B,, should remain continuous. In accordance 
with equation (7) it results in 


1 (22) 7 = ( 2) 1 (22) tie 
Fg, oy y=-0 oy y * cy y=d-0 Mz. oy y=d+0 


(=) (=) ; ( 4 ( ") (20) 
0x | y--0 Ox | y- +0 OX } y=-a-0 CX | y-a+0 


5. Solution of the differential equations 

Let us begin by solving the problem in the stator area. Considering 
equation (4) and boundary condition (18), we shall try the following ex- 
pression as a solution of the differential equation (10): 


2 
Ea = (y ‘i 3 dy) i Jw po by imax > V>Yy> dy. 
By, 7 om 


E= e-irat| A cosh 
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The solution of (11) in the air-gap area will become accordingly: 
E = e~?*'7( Bsinh(zy/r)+ C cosh(zy/7)]; d>y>0. (22) 


(Considering the boundary condition (18), we shall assume the solution of 
13) in the rotor-fixed system of coordinates to be 
E’ = e-i7@'"D cosh{n(y’ —d—d,)]; d+d,>y'>d. (23) 


r 


In the stator-fixed system it should read: 
; D ” 
E = e~i72l"— cosh[n(y—d—d,)]; d+d,>y>d. (24) 


The coefficient *%, which is normally a complex quantity, is to be computed 
as follows: the last expression, together with (13) and (17), yields 


-9 
(7, 7)? 3 Hy? 





[(L/o)(z/r)?+ 2rd,] — — = U(n/r)?, (25) 
eae 
52 re 2 SWLg Hy, F ‘ 
and hence ne — Pee (p/7)24 - __ J8OH6 Hr, _ 26 
. fly 1+ 2rd,o Ualr 7) (38) 


The solutions (21), (22), (24) and the boundary conditions (19), (20) pro- 
duce the following four equations for the calculation of the four integration 
constants A, B, C, D 


l 
A sinh!" d dy) (Hz, / In,,)\ = B, (27) 
V (Har, Hy,) \7 | 
od B cosh(zd/r)+ C sinh(zd/r)] = — D sinh(fd,), (28) 
T j Hz, 8 
A cosh di (May H)| —Jwpg by, (7 7)*Unax = C, (29) 
, D . 
B sinh(zd/r)+ C cosh(ad/7) = — cosh(nd,). (30) 
s 
Hence A = jorpg by, (7/7) tmax As (31) 


where A is a complex quantity depending on the slip and the machine’s 
dimensions; its full expression is given in the Appendix. 

From the equations (21), (31) one obtains the electric force in the stator 
teeth zone: 


E i Wo ML, Ya (7 77)*0 'max cosh” (y+ ; dy), He, Hn) - ; 


0>y>—dy, (32) 
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its mean value with regard to y being 


Bn y. es a E dy 
dy . 
—ads 


jralt 
€ JMLo- 


~ 





V1 


a sinh{(7dy/7),/(u,,/M@y,)} , 
(7 7)* max f d i — Sl. (33) 
1 (776 st T)y, (ary By) 
6. Stator’s overhang 
In all the above calculations we have completely overlooked the influence 
of the magnetic field on the parts of the conductors connecting the corre- 


p=oo 


Bye 
Pe ili 





Fic. 3a. Longitudinal cut of the machine with overhang. 


sponding stator slots. It is as though one were to neglect the overhang 
leakage in the classical theory. 

Fig. 3a shows a longitudinal cut of the induction machine. In order to 
represent the overhang on both sides of the stator, the rotor has been drawn 
in its natural length, while the stator has been lengthened up to 4J,, where 
l, is equal to one stator-turn’s mean length. Because of the great distance 
between the stator’s overhang and the rotor’s end rings, practically no 
magnetic flux links them. 

The stator overhang conductors are normally bent so as to recline on the 
stator’s yoke. This is taken into consideration by continuing the yoke, 
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behind the stator conductors, up to the length }/,. Fig. 3a shows also the 
relative permeabilities of the various regions, while Fig. 3b gives a separate 
view of } urt of t he dey el yped overhang alone. 


We are choosing a system of coordinates similar to that used in the 














Fic. 36. Part of developed stator overhang. 


previous calculations. The current density is to be assumed identical with 
that of (4). The differential equation of the intensity of the electric field E, 
ithe overhang zone is given by (10) with the sole correction p,, = p,, = 1; 
ence 

ek, @k, . , 
g . ae : JMpPo?. (34) 
0x* oy” 
the differential equation of EZ, in the air-gap (which now extends from the 
tator to the machine’s shaft, or, in the developed model shown in Fig. 3, 
wy=oo)l1s 


2B, @E e 
ot hel oe Oi (35) 


Ox? oy” 
The boundary conditions to be fulfilled by the solutions of (34), (35) are: 


i) Because of the perfect permeability of the stator’s yoke the tangential 
ignetic intensity H, should vanish along the plane y —d,,. In accor- 


lance with (7), this results in 


=") on @, (36) 
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(ii) The tangential magnetic field intensity H, and the normal induction 
B,, should both remain continuous on the boundary plane y = 0. In 
accordance with equation (7) one obtains consequently 


aD \ ap an an 

( =) _ (‘ =) ; ( =| = ( 4 (37) 
rose a ’ oy —— oe ° vi 
Ox y 0 ox y=+0 oy y=-0 oy y=+0 


Considering the first boundary condition (36), we shall write the solution 
of (34) in the form 


E, = e~I*#I"{K cosh[z(y+dy)/7]—jwpo(t/7)%maxti 9 > y > —dy; (38) 


maxs 


the solution of (35) will thus be 


f= ig Pew; yg > 0. (39) 


ee 
The last two solutions, together with boundary conditions (37), yield 
K sinh(zd,/7) = —L (40) 


and K cosh(zd4/7)—jwpo(t/7)*tmax = L- (41) 


By addition of (40) and (41), one obtains 
K = e-"4e'Thay19(7/77)*imax- (42) 


Inserting this value for K in equation (38) yields the intensity of the 
electric field in the stator overhang, 


f 
max 


E,, = e972" jwpyg(t/77)*imaxie- 74" cosh[m(y+d,)/7]—1}; 


0>y>—dy. (43) 
Its mean value with regard to y is 


7 7 ivgibess a _, _ sinh(zd,/7) 
+ = « E. dy — e-It: "wy (T i. e- a_i. etl! oF. 
d i 
st ; 70 4/T 


(44) 


7. Calculations of the stator current 

We assume the stator winding to be of the three-phase symmetrical type 
consisting of p pole pairs, q slots per pole and phase, and w conductors in 
each slot. All the conductors are connected in series. The e.m.f. induced 
in the conductors of the stator teeth region is added to the e.m.f. induced 
in them while running in the overhang zone. At the ends of each stator 


phase an e.m.f. E,,, is induced, where 
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_ sin(z 6) 


Eun = 2*pu 


sin(z 6q) UE max +1, (#a)max} ; (45) 
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here 1, is the length of the overhang on both sides of the stator (see Fig. 3). 


Inserting in the last equations the values of E,,,,, and (E£,)max a8 given by 


93) (44). we obtain 


sin(7/6) . 
7 94 9 olan ie 
Ey a*pw- Wo (7 7) bmax ” 


sin(7 6g)” 





sinh d tv \P Hy») | 
aa Lol FY y— TA yIT gy rd. 
= d 1 | +0, | ae a 1 . (46) 
| 71 dn (H,, My,) TO 4 | T 


Considering equations (3), (4), from the last expression an explicit relation- 
ship between the phase e.m.f. and the stator current J is obtained, 


9 


E Tjwpi, w . ae 
ph Jorg P 7d, 27 sin?(7/6q) 
sinh ” dy/(u,, Hy.) , | 
P 112 T aeeae rdyirSinh(7d.y/7) 1|\. (47) 
7 7 /T 
— dn /(Hz,/My,) 


The outside voltage U,,,, applied to each phase opposes this e.m.f. and should, 
in addition, cover the ohmic voltage drop / Ry, where Ry (Q) is the ohmic 
resistance of each phase. In accordance with equation (II) (see Appendix), 
and developing into a series the last expression of (47), we obtain 





17 pr ss 3jwpL,y wp 1. : +3(Ha,/Mr, dyn tanh(nd,) 


| © © 2a sin?(7/6q) (nd T)+(ady/th,,)+ (M7 u.,,7)tanh(nd,) 


+1,(] —8ndy/t+4(ady/7)?+ )| | (48) 


The last equation expresses the stator current in magnitude and phase 
as a function of the slip s. 

It is therefore an analytic counterpart of the well-known circle diagram. 
Fig. 4 compares the curve drawn according to (48) with the circle diagram 
of the same machine.+ It is significant that, whereas for small values of s 
both curves agree fairly well, the current given by (48), for s > 1, exceeds 
substantially (by about 40 per cent.) that given by the circle diagram: while 
at the same time the phase angle given by (48) is smaller than that shown 
by the circle diagram. The above differences arise because of the pro- 
nounced skin-effect, for higher values of s, which is neglected by the classical 
circle diagram. The rotor frequency is sw, and in the vicinity of s = 1 it 


The } s data and tl ircle diagram have been extracted from Say, Performance 
De 1.¢ Machines, pp. 392. 397 (Pitman, 1948). 
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approaches that of the stator. Due to the skin-effect, the current density 
at the outer part of the rotor is greater than (often several times) that in 
its deeper parts. Thereby the rotor effective resistance is increased and the 
stray magnetic flux decreased. The overall consequences are an increase 
of the rotor current and power factor. 


VU 





L 








$=0 


Fic. 4. Current-vector diagram (comparison with the classical circle). 


8. Case of no load 


When the machine runs at synchronous speed, i.e. zero slip, equations 
(48) and (II) (see Appendix) yield the no-load stator current 


, 3jwpyw*p | - ~ 
i, = 0, js | +1 [1—3(ad,/7)+4(ady/7)]}. 
0 on 5 sin2(7 6q) \d+ (dy by) + (de Hy.) | # | 5 (770 st ) 4 (776 st T) I 


(49) 


Ir/a 


R, in the denominator has here been omitted since it is very small compared 
with the other imaginary part. It is evident from the last equation that /, 
lags U,, by 90°, and thus no power is being transferred, a result obvious 
at no-load. 

[t should be borne in mind that J,, according to equation (49), is but 
the magnetizing current. The active component of the no-load current, 
which supplies the iron and other losses, must be calculated separately. 
Neither ,, nor «,, appears in equation (49) since there is no current at 
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all flowing in the rotor, and only a negligible one in the stator. The entire 
magnetic flux passes through the teeth zones in the y-direction and it is 


practically closed via the yokes exclusively. 


9, Case of infinite slip (s 1) 
The stator current is obtained from equations (48) and (III) (see 
Appendix) 
7 ) 3]WLy wp ( e on 
lo . oT Se sin?(q 6q) (ola 7d a" l| I 5 (md 7) + 4 (7d y 7)?]} 


(50) 


Neither the rotor’s permeability p,,, in the x-direction, nor p,, in the 
y-direction, appears in the last equation. At infinite slip the rotor frequency 
is also infinite. All of the rotor current flows through the outer layer 
exclusively. The magnetic flux, therefore, does not pass into the rotor 


teeth region. 


10. The torque developed by the machine 

This can be calculated in two ways: 

|. From equation (48) one can obtain the active stator current, and hence 
the power supplied to the stator; this after deduction of the stator J? R losses 
yields the rotor power and the torque follows. 

2. Equations (32), (7) give H, and E on the inner surface of the stator 
for y = 0) so that we can construct the complex Poynting vector on this 
surface. The real part of its total flux constitutes the power transferred to 
the rotor and the torque follows. 


Fig. 5 shows the curve of the torque (4) as a function of the slip s com- 
pared with that resulting from the circle diagram. It is worth while noting 
the difference of about 25 per cent. at standstill. For normal running 
conditions both curves yield practically the same results. 


11. Higher space harmonics 

Equation (3) gives only the fundamental harmonic of the stator’s linear 
current density which consists actually of harmonies 1, 5, 7, 11,..... To each 
harmonic a solution of the differential equations (10), (11), (13) corresponds, 
similar to that given by (21). We obtain therefore 1, 5, 7, 9,... space har- 
monies of E and H and consequently of the stator current. It can be shown 
(5) that these harmonics contribute only a few tenths of 1 per cent. to the 
stator current, and should therefore not be included. 


As to the torque developed by the higher field harmonics, it could be easily 
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calculated by constructing the complex Poynting vector of the relevant 
harmonics of F and H (4) as outlined in the previous section. 
All the above considerations have been based on a crude machine model 
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——— BY CIRCLE DIAGRAM METHOD 


Fic. 5. Torque/slip characteristic. (Comparison with the classical method.) 


which, while sacrificing the finer details of construction, provides a suitable 
basis for the mathematical analysis. It is clear, therefore, that the parasitic 
torques due to the higher field harmonics can be computed, while those 
caused by the intrinsic construction (dentated stator and rotor surfaces) 
are not revealed at all. 


APPENDIX 

Calculation of \ and related expressions 
(a) The quantity A, which is defined by (31) and calculated from (27), (28), (29); 

(30), is given by 
en tanh(7zd/7) +-("7/,, 7) tanh(’d,) ; (I) 
(Mx My,) 2 sinh{ (ad,/7), (Mr,/Py,)} + cosh (77d,4/7)4/(y,/y,)} tanh(ad/r) 4 

+ (7), 7)tanh(nd,)[(w,, Ly,) } sinh{(zd,, T)a/ (er, /y,)} tanh(ad/r)+ 
+ cosh{ (ard 54/7) (fez, /Py 
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Using the series coth a . 2 * 
_— he ee 
and the fact that in normal machines 
d, d; d,, > d; 4(,/My, (ard,,/7)? <1, 
ne obtains from (I) after neglecting the smaller terms, and inserting the very small 


argument 7d/r for tanh(zd/r), 


1) sinh 7754/7) a) (he, / pe 7rd 54 1+-4(prz,/z, d,, tanh(7d,) 
© (md g4/7)/ (pee, |p fey, T (wd /7) + (ad y4/pry, T) + (7i7/mz, 7) tanh(fd,) ° 


i 





(II) 
b) At no-load the machine runs at synchronous speed and the slip s is accordingly 


zero. Equation (26) takes the form 





n (77/7 )a/(ferg/Pys) 
The argument iid, = (7d,/T),/(j1z,/Hy,) IS Very small and can with reliable accuracy 
replace the tanh(7d,). 
As d(u,,/u,,)(7d54/7)(7rd,/7 1, expression (II) now takes the form 
sinh{ (ard,;/T),/(u Ly.)§ d, AL ° 
1—d VAP sve —_—___fn____.. with e = 0, (IIT) 
(1754/7) 4/(Hry/Pys) d+ (d5t/by,) + (4,/My,) 


) In the case of infinite slip (s > 0), % as given by (26) tends to infinity and 
expression IT) should read 

sinh{ (ad,,;/7r)4/(uw,,/py, )} “ : = 

1—A VAP rs Ps (Ly, /py, (7d,,/7)? with s 00. (IV) 


(ard ,4/T)4/(u by) 


The author is greatly indebted to Professor F. Ollendorff of the Hebrew 
Institute of Technology for his most valuable help and encouragement and 
to Professor R. Riidenberg of the Harvard University for the important 
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STATIONARY ROCKET TRAJECTORIES 
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SUMMARY 

The conditions to be satisfied by a rocket trajectory joining two specified terminals 
in space, if it is to be stationary relative to the fuel expenditure, are calculated, 
The rocket is supposed to be moving in a general gravitational field and to be 
subject to aerodynamic forces which are known functions of its velocity. The general 
theory is illustrated by application to the problem of the optimum programming 
of fuel expenditure to achieve a given height with a rocket ascending vertically 
from the earth’s surface, subject to uniform gravity and air resistance proportional 
to the square of its velocity. 


1. Introduction 
IN a previous paper (1) the fundamental conditions which must be satisfied 
by a rocket trajectory of least fuel expenditure joining two terminal points 
in vacuo have been calculated. Allowance was made for the existence of a 
general gravitational field, but aerodynamic forces which are functions of 
the rocket velocity were excluded. Tsien and Evans (2) and Hibbs (3) have 
studied particular cases of this problem of the optimum motion of a rocket 
subject to dissipative forces. It is the object of the present paper to extend 
the in vacuo theory in such a manner that it will deal with such cases. 
It will be assumed that the two terminals of the rocket’s path are given 
points in space and that the rocket’s velocity is known at these points and 
that its mass is known at the arrival terminal. In the author’s experience 
the majority of problems in this field may be so arranged that these assump- 
tions are valid. Forexample, Hibbs’s problem concerning the programming 
of fuel expenditure to achieve maximum range in horizontal flight over a 
flat earth may be posed in a form where, the range and the all-burnt mass 
being given, it is required to calculate the mode of fuel expenditure which 
will require a minimum initial mass. A number of given ranges might be 
selected and the corresponding optimum solutions calculated. The minimum 
fuel necessary to achieve any given range could then be calculated by inter- 
polation in these results. If, in any particular case, such a procedure proves 
to be inconvenient, there should be no difficulty in modifying the theory 
of this paper to suit alternative boundary conditions. 


2. General theory 


x; (t = 1, 2, 3) are Cartesian coordinates of the rocket with respect to 
fixed axes at time t. M is the mass of the rocket at this instant, c is its jet 
velocity, and /; (i = 1, 2, 3) are the direction-cosines of the thrust. The 


[Quart. Journ. Mech. and Applied Math., Vol. VII, Pt. 4 (1954)] 
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rocket will be supposed subject to two forces in addition to the motor 
thrust, (i) gravitational and (ii) aerodynamic. The gravitational force will 
depend only on the rocket’s mass, its position in space, and ¢ (allowing 
for a variation of the field due to the motions of planets, etc.). The aero- 
dynamic force will be assumed calculable when the rocket’s velocity and 
its position in space are both known. In practice this force will also be 
dependent upon the orientation of the rocket relative to its line of motion 
and upon the setting of various control surfaces, but we shall suppose that 
a particular orientation is automatically maintained during flight and that 
the setting of the control surfaces does not vary appreciably, being limited 
to the adjustments necessary to maintain the orientation. If this ideal 
situation is deemed to depart too radically from the actual conditions in 
any problem, the functions representing the aerodynamic forces must be 
allowed to be dependent upon additional variables and the theory below 
modified accordingly. No intrinsic difficulty will be found to present itself 
when we strive to achieve such generality, but the expressions involved 
in the statement of the theory become rather unwieldy. We shall there- 
fore take the external forces, to which the rocket is subject, to be func- 
tions of t, v;, #; (2 |, 2, 3), and M only (dots will denote differentiations 
with respect to t). The equations of motion are thent 


aM x : 
Mx; lic- nt F(t, x;,,%;, M), (1) 
where, in this equation, and in those which follow, j is to range over the 


values 1, 2, 3. Putting w clogM, —F,/M = f,(t,2;,%;,,w), we can 
write equations (1) in the form 


dw 7 , 

t dt “ this (2) 
3 

and then, since > 1, we obtain 
p=] 
dw lj 03 » ot 
(x,+f,)*}, (3 
dt A 4 , ty 


the positive root sign being appropriate, since w increases with t. The track 
being specified by equations x; = 2,(t), equation (3) is a first-order differen- 
tial equation for w. Our problem is to find functions 2,(t), satisfying the 
boundary conditions at the terminals, such that the increment in w over 
the trajectory is stationary. 


+ In equation (1), and in all later work, F';(t, x;, &;, M) is understood to denote a function 
Ve, Le, A 7 3 


the variables ¢, x,, 22, 23, €), Xo, Xs, 
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As explained in (1), the equations determining the stationary track 
become indeterminate when the equations of a trajectory of null thrust are 
substituted. From experience it is known that portions of a stationary track 
are often null-thrust ares. To avoid the exclusion of such a possibility, we 
shall modify equation (3) to take the form 

dw Tes ) 

dina Se Y (a.+f,)?+ 2 

dt PA fi) el (4) 
taking « = 0 when conditions on a stationary ‘powered arc’ are to be 
investigated and letting « > 0 when conditions on a stationary null-thrust 
are are to be discussed. 


We shall write 


/ : oe P\9 9 ° ee a 
LP) (x; +f;)?4 | = hit, x;, %;, %,, w), (5) 
1 ’ 
so that = nk (6) 
dt 
We now introduce a parameter A so that the equations 
xj= 2,0), t= tr) (1 


define a rocket trajectory and A = , at the departure terminal and \ = ), 
at the arrival terminal. The boundary conditions to be satisfied at the 
terminals are 


XL; = Xyq, L; = Ujp (at departure), (8) 
XL, = Ly, Lb, = Use (on arrival). 


We shall attach a subscript 0 to any symbol whose value is to be taken at 
the departure terminal and a subscript * to any symbol whose value is to 
be taken at the arrival terminal ; t, and t,, will be selected so as to minimize 
the fuel expenditure. If to, t, are given quantities not open to choice, or 
if the time of transit is fixed, some simplification of the general theory 
results. 

Employing primes to denote differentiations with respect to A, we have 


t= at, &; = (t’'xs—t"a')/t'3, (9) 

and equation (6) is equivalent to 
w’ = Hit, t,t", 24,24, 24, w), (10) 
where th = #. (11) 


Previous work indicates that optimal rocket trajectories involve im- 
pulsive thrusts at one or more points. In order to investigate the con- 
ditions to be satisfied in the neighbourhood of such a point, we shall allow 
the trajectory specified by equations (7) to possess one such point at A = «. 
At this point w is discontinuous by an amount equal to the magnitude of 
the vector velocity increment (1). If we employ a subscript — to indicate 
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that a quantity is to be evaluated at an instant immediately prior to the 
impulse, a subscript -+ to indicate evaluation at an instant immediately sub- 
sequent to the impulse, andif W denotes theincrement in w over the impulse, 


W= | s (é,,—2,_)*}. (12) 

Suppose the equations 
a ;(A)+-62,(A), t = t(A)+8t(A), w = w(A)+d6w(A), (13) 
specify a neighbouring motion to that given by the equations (7), satisfying 
the boundary conditions (8), and suppose that the fuel expenditure is 
stationary with respect to this variation, i.e. given that dw, = 0, then 


also Sw, = 0. dw is discontinuous at A x by an amount 5W, where, by 
differentiating equation (12), 
3. . . ae ed 
+a (14) 
( ~ a= ———--—— , 
| ¥ (a, _ 4 y2\ 
i ae y) | 


In view of equations (9), we have 


: t’dx.— a dt' " 
ox - — (15) 


J t’2 


> (x,t —ay_ t!, fr? (t, 8a; v; Ot’, )—t'? (t_ dx;_—aj_ dt_)} 
dv i a cg eee ee 
19 4ra | S r , J ’ o\4 
ral 1 2, (Xj t_—a; t.)*} 
_— (16) 


Substituting from equations (13) into equation (10) and neglecting 


second-order terms as always, we obtain 


oH,  oH,, cH, > (eH, oH, OH, .\ , OH 
ow ot 4 — ot - ot S —— 0%;-+ —,; 02%,-++ —, 02;)}-+-— 6 
ot ot ot L~ \0x; OX; OX; ow 
a j j j 
(17) 
This is a first-order linear equation for 5w, which we integrate over the 
interval (Ay, x), taking as the initial condition dw = dw, when A = Ag, and 
obtain 


eno { 2a 


x 


bu, | exp| : * 0H an)|— a pa St’ +O at's 5t”- 
¢ 


J cw \ ¢ ot’ 
L 


Ao A 


— 
a (= mrt "re cH 5 bx) aa}. (18) 
P. ox; * ea” 


— 


_ 
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An er by parts shows that 


[ exp( "oH ale OH sy ay 
ot 


|) éw 


. eH ,.\eH . ie "oH .\oH 
= exp| — f oa) ay =| ines — | — an)? ‘a dh 
lo A 


Xo ; Ao 0 Ao 


—— se 


om [hal 





- * oH eHeH d (oH 
ll | exp( — J ew a)[# ‘a ala) | ——_ = 
Ao : Ao 


and, similarly, two integrations by parts give 


{ esol — | IE ra st” dA 


0 


= exp(— | _~|j— H a) OH dt’ + es oS es . = stl — 
; | at ew et dA\ et”) | _ 


Ao 


| 
a 


° 7 oH an)| oH 90H d eH = 
7 fex(- [ ew oo) et” ew drd\ at’ 


_ oH dad oH\ 2 oH stdr. (20) 
ot” drA\ew] | dd2\ at” 


_ {eH Se’ oH oll d pla 


| ot” ew et” 


There are six more results similar to equations (19) and (20), decid from 
{ 


them by replacing ¢ by the 2. 
Equation (18) can accordingly be written 


"oH ..\r. ; ee 
éw_ = exp| | . aa) su. A 9 dty— Dy dty — > (Bip Sx jo + E59 dX}9)|+ 

J éw j=1 

Ao 


3 
+A_8t_+D_8t'_ + > (B;_8x;_+ E;_8a;_)+ 
j=1 


H an)(7 66 | ¥ X,8x,) da, (21) 
j=1 


whe 
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where 


CH cH cH Ser) 
s at’ | dw at” da\. 


OH cH cH d(ceH 
st dd cx, Ow ex; dr\K x; 
teas f E, on 


at” da\at’ 


at | cw ot 


T 0H oH eH (—) oH ner 


\OW | 


,0H d(eH\ eH d(éeH\ , d® (eH 
“ew da\ at’) at” dA\éw] * dd2\ at” 
CH oH eH CH\?*eH d{(eH 

H| ) ex, dd 


, oe >. 
Ox OW CX Ow OX; 
j / J 


= nl) = (Fa) 4 aaa) (22) 


dw dr\éex' Ox; dX\\ew] © dd? ex; 


By integrating equation (17) for dw over the interval (a,A,), we obtain 


similarly 


3 
iw, = A, dt,+D, dst’, + > (B;,8x,;,+H;, 82;,)— 
j=1 


de 
20 exp| - oH - dr)| A 


x 


+ dt, +D Witt S (Bye dt +Ey. Sarj)|— 





| feso iE \}( 7 3t + ¥ X,dx,) da, (23) 
j=1 
| 
nee dw, 0 by hypothesis. 
Since, however, the discontinuity in dw at A a is of amount 6W, we 
Sw,—dw SW, (24) 


ind hence, substituting from equations (21) and (23) into this latter 
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equation and solving for dw», we obtain 
Oy’ 3 iw] ’ 
SW) = Ay dty+ Dy Sty + 2 (Bio 82x jg + Ej9 8X59) + 
{= 


( eH 


ew 





+exp| — 


© 


Ao 


da)[A,.8t,—A_8t_+D, St’. —D_8t’. + 
3 
+ > (B,,82;,—B,_82,_+E,,3%;, —E,_ 8a )—3W|— 


gu} 


oH oe pe . 7 oe 
—exp| _— E — a)[A. dt, Dy Bt +X (Bix Ditjx Ei 4 32ix)|— 
F. 


ow 
Xo 
7 H 
~ | exp| _ l a an) T n+ X ; dx;) da. (25) 
Xo } a 


In order that the varied path shall satisfy the boundary conditions, we 
require that 82;) = 62;, = 0. Also, 5% = d%;, = 0 and hence, from 
equations (15), we may conclude that 


$2/, = ( a , Oty = (7 a ; (26) 
t 0 t * 


dt, and dt, will also be zero if the times of departure and arrival are fixed, 
but we shall suppose these times not to be specified. Since the 2; and t are 
everywhere continuous, d2;, = 52; Az;, dt, = d5t_ = At. Otherwise 
the varied trajectory is arbitrary. 

Equation (25) can therefore be expressed in the form 


3 ’ 
. a i 
dW (p— > B72) Bto-+ Ay Stat 





J 
j7=1 
Hexp| io (= aa)[(A,—A_)At+) S (B,. —B,_) Az;+ 
Xo : 
3 
Dt. —D-3F. > (E;, 8x; E; Ou; _) aW| 
j=1 
* oH S. : 
> j— Le ot, | — 
exp | = an)|(z 2, Bi), 1, dt, 
Ao a 
*. : . 
- | exp| | = an}(T 31 + ¥ X,8x,) dA, 27) 
ow j=i J ) 


Xo Ao 
dW being given by equation (16). 
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If the path specified by equation (7) is a stationary one, 5w is zero for all 
- . = vor ~ OQ,’ Oy . ¢.,/ / a, a, 
values of the arbitrary quantities dtp, dt,., dt, ty, Ax;, At, a; , daj;_, dt, dt_, 
and for all forms of the arbitrary functions 52;, dt. The coefficients of these 
quantities in equation (27) are accordingly all zero. Thus, 


A, = A, = E,.—P. = D_—R_ = E,,—B,, = D,—B, 


J 


= j+ = 
X, T = 0, (28) 
where — —" 
’ x;,t_—az;_t P 
J u > 
(P= a =e 
> (2; t v5 t "| 
1 
. _— : 
> 2; (2; t uj t',) 
t’2 R j = a ssiaisss 
} ne / , , , »\? 
| 2 (24! vj t.)*} 
j=1 
3. , , , , , 
> xj. (x; t x; t.) 
tt R, ~ Paes t (29) 
> (x54. ¢_—aj_t',)?} 


_ 


, 


Lt = © (30) 


3 
D—NE 
in * 


j=1 


In addition 


it both end points. 

This latter condition is easily shown to be satisfied at all points of a 
trajectory, for, by differentiating equation (11) and employing equations 
9), we obtain 


cH Ss ah di, 3. x! ah 
D S (aes } Oe, (31) 
Ct hone Cx; Ot ft" On; 
| @H th et, él 
E , as OF aed Be, (32) 
, O2 OX; OX; t CX; 


It is now clear that equation (30) is an identity. 
Further differentiation of equation (11) and appeal to equations (9) leads 
to the following results: 


cH Ch cx: ,oh 0« ch t” oh 
. “+t — 4 : ~ 


dors ; 5 ; (33) 
Cx OL; Ox; CX; OX; Cx; ts Cx; 
cH h : Wee OX; , y oh a] 
ct hol ex; ot 0%; ot 


. Oh , 3t "2; 
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And now, by reference to equations (22), we may deduce that 


3 ae oe 
hh... oh oh oh d {ch 
BS on |. @ Lz. +4; fae ie oie ata 36 
p3 "aa, * 8, o * 0x; Ow ? sae} a. 
J 


BR. = ch P Oh oh ; waz) 


+ - 37 
7 02; Ow 0%; dt “ 


OX; 


Substituting for h from equation (5) into equations (31), (32), (36), and 
(37) shows finally that 


'D = —r x4 (é+f), vB; = (@+f)r, (38) 
[So oe 
where r = |, tf) +e?) , and 
3 
_— I ve of; | — 
Aas > | far - > af) +¢] 
~e > 4 (+f) > (iy tf) Ze > TEA) (39) 
3 , 
_ t. (Sie, dk ee \ d ' 
B,; = > |¢ ff Nex +r ok a; f\|- —F(rlés+f))}. (40) 


Equipped with the above expressions for A, B,, etc., we can now apply 
the conditions (28) in any particular case. We will now summarize these 
conditions. 

(i) If the times of arrival and departure are free to be chosen so as to 
minimize the fuel expenditure, then A must be zero at both end points. 

(ii) Across an impulse, the quantities A, B; must be continuous. 

(iii) On either side of an impulse Z; = P;. Together with equations 
(29) and (38), this implies that 

eth 054 Hj - (41) 


3 , 
\! | . - \2\” 
(S x; itt) +e} [Se -#,) | 


on either side of an impulse. If ¢ = 0, the left-hand sides of these equations 
are the /;. The right-hand sides are the direction-cosines of the impulsive 
thrust. These conditions may therefore be expressed by saying that the 
direction of thrust must be maintained steady from one side of an impulse 
to the other, i.e. the direction of thrust is continuous through an impulse. 
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The case where an impulse operates at an end point of an arc of null thrust 
will be dealt with later (section 3). 

The conditions D R_, D, = R,, will be found to be linear combina- 
tions of the conditions (41). 

(iv) In addition, the coordinates x; must be continuous across an impulse. 

Along a stationary arc, it is necessary that 7 = 0, X;= 0. These 
equations may be written in the form 


eH eH, dA 


T’ ade ne ee | 
ot ow dxr 
, H oH dB, 
X | ae ee | 42 
; 6x; Ow B; dX (42) 


Only three are independent (4). Selecting the latter three and taking A = ft, 
we have, along a stationary arc, 
h . ak 
ee a. (43) 
OX; OW dt 
whence, substituting for h from equation (5), we obtain the equations 


3 


ey * tf (fe, S| _ 4B; = 44 
ps Tk ess Biol dt 0, a 


the B; being given by equations (40). With « = 0, these fourth-order 
equations for the x,(¢) and equation (4) are sufficient to determine stationary 
powered trajectories, apart from the values of certain constants of integra- 
tion which must be selected to suit the boundary and other conditions 
summarized above. When « + 0, these equations possess another type of 
solution, each of which approaches a null-thrust arc as « > 0, Such solutions 
will be investigated in the next section. 


3. Null-thrust arcs 


Along an arc of null thrust, 


x,+f 0, w = constant, (45) 


J vJ 


and, with « = 0, the equations (44) become indeterminate. Such an are 
is clearly minimal, however, since any neighbouring trajectory can only 
be followed by consuming fuel. This type of are may therefore be utilized, 
together with powered arcs governed by equations (44), in the construction 
of a stationary trajectory, and it is necessary to calculate the form the 
conditions (28) take for such an are. To do this we shall find the limiting 
form of these conditions as « > 0. 
5092.28 K k 
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If c is small, but not zero, equations (4) and (44) may be shown to possess 
solutions of the form 


u; = X;+;, w= We, (46) 
where the X; and W are solutions of equations (45) (i.e. W is constant) 
and we shall neglect terms O(e?). Such a solution approaches a null-thrust 
are as e« > 0 and is therefore appropriate for study in this connexion. To 
show that such a solution exists and to calculate the equations which must 
be satisfied by the functions €;, we first observe, if we substitute from 
equations (46), that 


t;+f; = Xj +E; +hi(t, X,+e€;, X;+é,, W+ef) 
= X,+f;(t, X;,X;, W)+ 
= é é 
sin S (eos Laos 
= €p;, (47) 
7} ee : Sf; we of; 
where p; = €;4 pA ( X bat ax E+ oy aW im (48) 


Hence, if we substitute from equations (46) into equations (39) and (40) and 
retain terms of zero order in ¢ only, we find that the limiting forms of the 
quantities A and B; are 


4=)> u(f)— S i, 2¢)- 
k=1 ~e 


fos j=1 
fe. A 
B, = 2, oloe + ap) —%; (50) 
3 a 
where u; = p;| > pg+ 1) . (51) 
K=1 
3 wl 
and hence p; = ull—> uz) . (52) 
e=1 


The limiting form of equations (44) are found to be 


3 Z 
— } k i‘. of, dB; = 53 
Zz (mae B; uy: nF = (53) 


k=1 
and of equation “a is 


3 
=f = | BY p24 1) - -(1 — > ui) ‘ (54) 
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the quantities B; appearing in equations (53) being those defined by 
equations (50). 


3 
Equations (52) show that, as } uj > 1, the p, tend to infinity and hence, 
(=1 


k 

in view of equations (48), the €; also tend to infinity, but the &; and &; 
3 

remain finite. As we approach a point where > uj = | therefore, the <; 


R71 

tend to infinity. Since the €; must be real, the u; must everywhere satisfy 
> uj ¥ 

[=1 

Equations (41) are to be satisfied on each side of an impulse and, ife + 0, 


i 


this implies that #;-+-f; must become infinite as we approach the impulse 
in such a manner that 


(4 4—2,_):(%Z_,—o_): (%g,—X3_). (55) 
But equations (47) imply that, over an arc of the type we are considering, 
(%, +1): (ote): (Fe+fg) = Py: Po: Ps = Uy: Ue: Us. (56) 


At an impulse on a null-thrust are it is therefore a requirement that 


> uz a Uy i Ug! Ug (4 4—2,_): (L_4.— Ze_) : (%g,—Z3_), (57) 


ie. the wu; are the direction-cosines of the direction of impulsive thrust. If 
two null-thrust ares meet at a point of impulsive thrust, this condition must 
be satisfied on either side of the point. 

The restrictions upon the values of the quantities A and B; at an impulse, 
deduced in the previous section, must also be satisfied at an impulse 
situated at an end point of a null-thrust are. On such an arc these quantities 
are calculable from equations (49) and (50). It is obvious that, if two null- 
thrust ares meet at a point at which no impulse is applied, they represent 
the same solution of equations (45) and need not be regarded as distinct. 
It may also be proved that the functions uw; on the one are are analytic 
continuations of these functions on the other are. We may therefore ignore 
this possibility ofa junction between two stationary null-thrust ares at which 
no impulse is applied. At a junction between a null-thrust and a powered 
arc an impulse may or may not be applied. If no impulse is applied at 
such a point, the quantities A, B;, x;, #; must be continuous there, but 


no other conditions have to be satisfied. 
4. A first integral 


In practice, it is frequently the case that the functions f, are not explicitly 
dependent upon t. If such is the case, a first integral of the equations (44), or 
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of their counterparts for null-thrust arcs, equations (53), is easily obtained, 
For H is not then explicitly a function of ¢t and hence the first of equations 

9 av Wr > 
(42) may be written oH dA ‘ » 
ew dd’ "7 


from which we deduce that 
iin Kexp| [ . aa), (59) 
| ow 


K being a constant of integration. 
The limiting null-thrust form of the first of equations (42) is found to be 


3 ‘ , 
of; of;\ dA 
Blatt. A Ht). == @, 60 
> (2 zt) dt oy) 
j=1 
from which we deduce that, if the /; are not explicitly dependent upon t, 
. 3 c 
. of; 
A = K exp! | > u, di dtl. (61) 
lJ \. ow) | 
j=1 
But A = 0 at either terminal and hence K = 0 on either arc leading 
from a terminal. Hence A is identically zero over both these arcs and, in 
view of its continuity, it is therefore zero everywhere. 
Accordingly the first integral is 


A = @, (62) 


5. Rectilinear motion 

The programming of the fuel expenditure to achieve maximum height 
with a rocket ascending vertically from the earth’s surface is a rectilinear 
problem of the type we have been discussing and has been considered in 
detail by Tsien and Evans (2). In any such problem the rectilinear motion 
may be assumed to take place along the z-axis and, omitting the subscripts 
which have served to distinguish between the three coordinates, force 
components in the directions of the three axes, etc., the conditions to be 
satisfied by a stationary track are as follows: 

The track will be composed of null-thrust and powered arcs, the latter 
satisfying the differential equation 


f , Fp dB _ 


7 ae — = 0, (63) 
~~ da Ow dt 
where B= 47 + F (64) 
0% Ow 


and we take the positive signs at the ambiguities in the case of motion in 
a positive direction and the negative signs otherwise. These equations are 
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the one-dimensional form of equations (44) and (40). Along a null-thrust 


are x = X(t), w = W, a function uw is specified by the equation 
c f of dB ~ 
(] : pm OF, coe B — — = QV, 6: 
‘ax “aw dt ae 
du . 
where B u 3 4 wT - % (66) 


i.e. the one-dimensional form of equations (53) and (50). Over such an are 
it is also necessary that |u| < 1. 
At any junction between two arcs the quantities A and B are continuous, 


A being given by the equations 


hws (/- #2) 27, (67) 
: \ Cx CW 
A u (f X fe) —y2X a + Xiu (68) 


on powered and null-thrust arcs respectively. 

An impulse may be applied at the end of any such arc, provided its 
direction is that of the immediately preceding or succeeding finite thrust, 
or, if the are is of the null-thrust variety, provided uw = +1 on this arc at 
the junction if the impulsive thrust is in a positive direction, and u = —1 
at this point if the impulsive thrust is in a negative direction. These are 
one-dimensional forms of conditions already obtained in section 3. 

If, as is generally the case, we wish to minimize the fuel expenditure with 
respect to the times of arrival and departure, then A = 0 at both these 
instants. Also, if f does not depend explicitly upon t, then A = 0 identically 
over the whole track and this equation is an alternative to equation (63) over 
a powered are and to equation (65) over a null-thrust arc. 

We shall conclude the paper by applying the above conditions to the 
case of a vertically ascending rocket subject to a uniform gravitational 
field and an air resistance proportional to the square of the velocity. Thus 
we take kei 

f= 9+— = 9+kz? exp(w/c). (69) 
M 
The initial velocity in the launching tower and at the height to be attained 
will be assumed to be zero. We shall minimize the fuel expenditure with 
respect to the time of transit between these points, so that A 0 every- 
where. 

Assuming the motion to be characterized by an initial climb under motor 

thrust from 2 — 0, followed by a subsequent null-thrust phase during which 











502 D. F. LAWDEN 


the rocket coasts to a standstill at the specified maximum height, the first 
part of the trajectory is governed by the equation 


A= 1—S— = 0 (70) 
(from equation (67)) and the equation of motion 
c dM : ka? - 
Ma TIT (71) 


Equation (70) shows that an initial impulse must be applied to the 
rocket to boost its velocity up to a value cV, where 
3s. “9g . 
y34/y? — —. (72) 
ke2 
This impulse could be provided by a booster unit which falls away after a 
few seconds of thrust. The subsequent motion of the rocket is calculated 
by eliminating M between equations (70) and (71) and integrating to give 


q 79 9 ry 2.1 S 
«x = 4(V?—s*)+2(V—s)+2log 
ce 2+ y 
"119 
I, — V—e-+- log! Of (V+2)\ (73) 
c | s(s+2) | 
where cs = #. The mass M during the climb is given by 
J 
ae. Met). (74) 
ke2 
Over the null-thrust are the equation of motion is 
- kxX2 
X+ 9-4 0, 75 
g M, (75) 
M,, being the all-burnt mass. Hence 
wX = cot(wt+¢), (76) 


where w? = kg/M, and ¢ is a constant of integration. We shall shift the 
origin of time so that ¢ = 0 at all-burnt. Then @ = wt+¢ increases from 
¢ to 4m as X decreases from its all-burnt value to zero. 
Since A = 0 over the null-thrust arc, equation (68) yields 
du ae . -- 
7a w(tan 6 — cot 6)u = u? cot?é. (77) 


The solution of this equation is 


( ke? ' sin 20 ‘ 
u : — (15) 
J shat — logsin?é-+ l 


U being a constant of integration. 
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B must be continuous across the junction of the powered and null-thrust 


ares. Over the powered arc, equations (64) and (67) show that 


B = (f—A)/# = fiz, (79) 
ad equations (66) and (68) similarly show that over a null-thrust are 
B = (uf—A)/X = uf /X. (80) 


If an impulse be applied at the junction and ¢_, X., denote the velocities 
before and after the impulse respectively, since wu = +1, the continuity 
of B requires that 

g+kz?_/M g+ kX2./M 


ees (81) 
But (X ,—a_) represents the increment in w over the impulse and hence 
M,/M exp/ = x, = } (82) 
Condition (81) is therefore equivalent to 
ie v)+vV(v—Ve"-) = 0, (83) 
where 7 z_andcV X ,. But equation (74) is satisfied at the junction 


on the powered arc and hence 
M_<¢ ; 
SF = v%v+1). (84) 
ke? 
Equation (83) may now be written 
y2(v 1)( V -v) nl wVi(v - Vel v) 0. (85) 
The impulse must be in the same sense as the thrust on the powered arc and 
hence V > v. Taking this inequality in conjunction with equation (85), it 
is easy to show that V = v is the only admissible solution. Thus there can 
be no impulse at the junction. 

[t now follows that #, and hence f, are continuous across the junction and 
the cont inuity of B (as given by equations (79) and (80)) therefore requires 
that w | at the junction. Thus when @ = 4, u 1. 

At the junction on the null-thrust are 

wX = cotd. (86) 

The velocity being continuous at this point, s = cot d/we at the end of the 

powered arc. Substituting in equation (74), we find that the all-burnt mass 
is related to ¢ by the equation 

V(M,,.g/ke*) tan*d — tan ¢. (87) 


Employing the condition u 1 when 6 = ¢ to determine U in equation 
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(78) and substituting for ,/(.M,.g/kc*) from equation (87), we deduce that 


2(1—z,)! z4(1—z)t 
ee aT o) : “ | ) - (88) 
23(2z)— 1) z—log z+(2—3z )(2z2,— 1)-1!+- log zy 





where z = sin?0, z) = sin*d. 

Equation (87) shows that tan¢ > 1 and hence z, > 3. As the rocket 
coasts from the junction to the peak of its trajectory, z increases from z, 
to 1. It will be found that w then decreases steadily from unity to zero 
and hence |u| < 1 over this arc, as is necessary if the trajectory is to be 
stationary. 

This mode of programming the fuel consumption therefore leads to a 
stationary trajectory. 

To make the above calculations yield results of practical importance, it 
is of course necessary to allow for the decrease of air density with height. 
This may be achieved with reasonable accuracy by taking the air-resistance 
term in equation (69) in the form kz? exp(—yx)/M and then repeating the 
above calculation. The reader interested in the results obtained, and the 
practical significance of these results in the design of high altitude rockets, 
should refer to the work of Tsien and Evans. 
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SUMMARY 
The double Fourier series expansion of the output from a half-wave linear rectifier 
with a two-frequency input is well known. The purpose of this paper is to provide 
tables of numerical values of the first several coefficients in that expansion as 
functions of the input amplitude ratio k. The tables provided are of sufficient 
accuracy to permit a number of additional coefficients of higher order to be deter- 


mined by known recursion relations. 


1. Introduction 
Iris well known that the output from a half-wave linear rectifier responding 
to a two-frequency input 

a(t) = P[cos(pt+-6,)+k cos(gt+6,)], 9< k <1, (1.1) 


can be expressed as a double Fourier series in the form 


x 


y(t) = $PAg(k)+P S* Asmn(k) C08(M.mnt+bsimn): (1.2) 


m,n=0 — 
where the asterisk indicates that the sum is extended over both signs if 
mn + 0 and over only the upper sign if mn = 0 with the single exception 
that the term occurring under the summation sign for m = n = 0 is not to 
be included in the sum. In the expansion (1.2) the angular frequencies and 


phase angles are given by 


w.», mp+nq, Psmn = mO,+n86, (m,n = 0,1, 2,...), (1.3) 
and the amplitudes are the functions of k defined as 
< ae 
A, (k) , (cos u+k cos v)cos mu cos nv dudv (ns, 00 <= ©,.1,3...3, 
I (1.4) 


+ Presented by title to the American Mathematical Society, 23-24 April 1954. 


[Quart. Journ. Mech. and Applied Math., Vol. VII, Pt. 4 (1954)] 
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where R is the region of the square 0 < u,v < winwhichcosu+kcosv > 0, 
For a detailed proof of these facts we refer to the original work of Bennett 
(1), (2) or to the more recent paper of Sternberg and Kaufman (3). The 
purpose of the present paper then is to provide numerical tables of the first 
(k) = A. ny(k) as defined by (1.4) for 
k = 0-02(0-02)1-0 which may then be used in conjunction with interpola- 
tion for direct applications or may be combined with the aid of recursion 
relations to obtain values of the higher-order functions (1.4). Since the 
functions A,,,,,(k) = A..»,(k) are zero for all values of k for m-+-n odd and 
greater than one, the functions thus tabulated are $A (k), Ayo(k), Ao, (k), 
A x(k), Ay,(k), Aga(k), Ago(k), Agi(k), Agel), Ay3(k), and Ay,(k). Each entry 


in the tables is believed to be accurate to somewhat better than one unit 


eleven non-zero functions A 


<“mn 


in the last place given and most of the entries are thought to be rounded 
correctly. 

The construction of the tables is described in section 2, some brief 
remarks about the use of the tables are given in section 3, and the tables 
are given thereafter. The functions A,,,,(k) = A.,,,,,(k) are the special case 
(h,k), briefly tabulated by Sternberg (4), 
in which the bias h is zero. The name Bennett functions has been proposed 


of the more general functions 4A.,,,, 
in (4) for these more general functions which, in view of the relation 
A,,,(k) = A,,,(0,4) just mentioned, thus also applies to the functions 
under consideration. 


2. Construction of the tables 


As the methods employed in constructing the tables were such that it 
was about as easy to prepare an eight-decimal table as to prepare one of 
lesser accuracy, it was decided at the outset to provide this relatively high 
accuracy even though direct applications at the present time may require 
less. One immediate advantage of the high accuracy maintained will be 
noted in section 3 in connexion with the use of the tables and recursion 
relations to obtain values of the higher-order functions 4,,,,,(k). 


m n(k) were evaluated 


To begin with, for 0:02 < k < 0-40, the functions A 
from the series 


4A oo(k) = (1/7)[1-+(1/4)k2-+ (1/64)k4+...], (2.1) 
A(t) = 1/2, Agi(k) = &/2, (2.2) 
A go(k) = (1/327)[2—(3/2)k?+ (9/32)k4+...], (2.3) 
A,,(k) = (k/7)[1—(1/8)k?—(1/64)k4—...], (2.4) 


Aoo(k) = (k?/32r)[(3/4)+(1/16)k2+ (9/512)k4+ ...], 
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TABLES OF BENNETT FUNCTIONS 
A,(k) = (1/157)[—2+ (15/2)k2—(225/32)k4+ ...], (2.6) 
Ag (k) = (k/37)[—1-+ (9/8)k?—(15/64)k4—...], (2.7) 
Ago(k) = (k?/m)[—(1/4)+ (1/16)k?+ (5/512)k4+-...], (2.8) 


A.3(k) = (k3/32)[—(1/8)—(3/128)k?—(9/1024)k4—...], (2.9) 
A o(k) = (k*/157)[(5/64) + (9/256)k?+ (75/4096)k*+...], (2.10) 


using from nine terms in (2.1) to five terms in (2.10). Each of these series 
converges for all k, 0 < k 1, and may be derived by direct integration 
of the inner integral in the corresponding double integral (1.4) followed 
by expansion of the integrand of the resulting simple integral in powers of 
keosv and termwise integration. For details of a typical derivation we refer 
to the derivation of the more general double-series expansions given by 


(h, k) 


Sternberg and Kaufman (3) for the functions of two variables A 


mn 


mentioned in the introduction. 


For 0-42 < k < 0-98 the functions A,,,,(k) were evaluated from the 
formulae 
4A yg(k) = (2/?)[2E(k)—(1—k*)K(k)], (2.11) 
A(k) 1/2, Ay, (k) k/2, (2.12) 
Ao(k) = (4/922)[(4— 2k?) E(k) —(1—k*) K(k)], (2.13) 
Ay (k) = (4/37?k)[(1 +k?) B(k)—(1—k*) K(k)], (2.14) 
Ayo(k) = (4/97?k?)[(—2+- 4k?) E(k) + (2—5k?+ 3k*) K(k)], (2.15) 
Ay(k) = (4/2252?)[(—38-+ 88k2— 48k) E(k)+ 


| (23—47k2+24k4)K(k)], (2.16) 

i(k) = (4/452?k)[(3— 13k?+ 8k*) E(k)+ (—3+7k?— 4k) K(k)], (2.17) 

Ago(k) = (4/1522k?)[(— 2+ 2k?—2k4) E(k) + (2—3k?+ i) K(k)], (2.18) 

A,,(k) 4/4522k3)[ (8 — 13k2+ 3k4) E(k) + (—8+17k*?—9k4)K(k)], (2.19) 
o(k) = (4/22522k*)[(— 48+ 88k?— 38k) E(k) + 

+-(48—112k?+ 79k4*— 15k*®)K(k)], (2.20) 


where K(k) and E(k) are respectively the complete elliptic integrals of first 
ind second kind with modulus k as tabulated by Fletcher (5). The formulae 
2.11) to (2.15) have been given by Bennett (1), while the remaining 
lormulae may be obtained readily from the former with the aid of the 
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recursion relations for the functions of two variables 4,,,,,(h,k) noted by 
Sternberg and Kaufman (3) and, according to Bennett (2), originally due 
to 8. O. Rice. These formulae are repeated below in section 3 with h = 9 
for convenience in connexion with the use of the tables. 
Finally, for k = 1-0 the functions were computed from the formula 
8(—1)” 


im in=—0,2,4,...) (5 
mn(1-0) [(m+-n)?—1][(m—n)?—1]r? idinitilieiaea 





bo 
Li] 


of Bennett (1). The actual arithmetic was performed on a thirteen-digit 
floating decimal I.B.M. Card Programmed Calculator and the results were 
rounded to eight decimal places. An examination of the computation 
indicates that each entry in the tables should be accurate to better than 
1x 10-® units and that a majority of the values are probably rounded 
correctly. This then completed the table proper. 

In an effort to ensure accuracy several checking operations were per- 
formed. To begin with, in an overlapping region surrounding k = 0-40, the 
tabular entries were computed to full accuracy on the thirteen-digit floating 
decimal I.B.M. calculator using both the series and the elliptic-integral 
formulae listed above. For 0-02 < k < 0-98 the entries were computed to 
a varying degree of accuracy using a separate eight-digit I.B.M. installation 
with both sets of formulae and the results were compared with the accurate 
entries. Finally, every fifth entry in each column of the tables was com- 
puted to full accuracy on a ten-digit desk calculator. As a last precaution 
the tables were differenced after typing and an examination of the differ- 
ences appears to substantiate the statements concerning accuracy made 
above and in the introduction. 

We note in passing that the values of A,,(0-02) and A,,(0-04) given in 
the table as zero to eight decimals are more accurately A,,(0-02) = 3x 10-" 
and A,,(0-04) = 4x 10~-® respectively. 


3. Applications of the tables 


For direct applications of the tables linear interpolation should be 
generally good to at least four decimal places in all of the functions A,,,,,(k) 
tabulated, with the accuracy gradually increasing to the full eight decimal 
places in the case of A,,(k) with k small. If greater accuracy is desired, 


mn 


interpolation with the Gregory-Newton formula using second differences 
may be employed and should yield results generally accurate to about six 
or more decimal places for k < 0-80. If still greater accuracy is desired, a 
five-point Lagrange formula is suggested and should give results accurate 
to about 2x 10-* units for k < 0-80. 
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For purposes of obtaining values of the higher-order functions A,,,,(k), 
the high accuracy of the tables permits them to be used to advantage in 
numerical evaluations employing the recursion identities 


(m—N+-3)A nig na = —(M4+N—3)A py _yn—1— 2MkA yy, 
(m+n+1)A,,, (m—n—3)A po ,—2(M—1)kKAp_yn—y, (3-1) 
(n+m-+1)A,, —(n—M—3) Am n-2—2(N—1)(1/k) Aga n— 

n—M+3)An in = —(n+mM—3)A ps n—1— 2M(1/kK) Amn: 


The identities (3.1) are equivalent to those of Sternberg and Kaufman (3) 
with the bias h = 0 and may be derived in the manner briefly described 
therein. Despite the advantages of the high accuracy of the tables just 
mentioned, the user is nevertheless cautioned to note that the factor 1/k in 
the last term of the third and fourth identities (3.1), together with the ever- 
present round-off errors in the eighth decimal place of the functions A,,,,,(k) 
tabulated, imposes a serious limitation on the accuracy of the values of the 
higher-order functions A,,,,(k) computed in the manner indicated above 
when k is small and n > m. 
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0:02 0°3153 4172 
"04 *3184 3722 
"06 *3185 9643 
08 +3188 1939 
‘10 3191 0616 
“12 | *3194 5084 
“14 *3195 7152 
16 =| = *3203 5035 
‘18 =| +3208 9346 
+20 | #3215 o102 

| : 
"22 *3221 7323 
‘24 ~=2+| ~=«*3229 1029 
+26 *3237 1245 
28 3245 7995 
-30 +3255 1308 

| 
32 02«| «3265 1216 
“34 *3275 7752 
+36 +3287 0952 
*38 +3299 0855 
“40 | °3311 7504 
"42 "3325 0946 
“44 3339 1228 
“46 +3353 8404 
“48 +3369 2531 
*50 *3385 3070 
"52 *3402 1857 
*54 | *3419 7253 


“60 *3476 7034 
62 3497 1818 
64 +3518 4196 
-66 *3540 4281 
68 } 3563 2195 
*70 *3586 5070 
| 
72 "3011 2053 
“74 *3636 4304 
76 +3062 4999 
78 =| = +3689 4335 
‘80 °3717 2531 
“82 *3745 9836 
84 3775 9532 
“86 | +3806 2943 
“88 *3837 9452 


4979 
4666 


*3904 


*3939 


“96 | °3975 7430 
98 | *4013 4409 
1-00 "4052 8473 
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‘03 
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“47 
48 
“49 
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Ago(k) 


o°2121 
*2119 
“2116 
‘2111 
+2106 


*2099 
*2090 
*2051 
*2070 
*2058 


"2045 
*2031 
"2015 
*1999 
“1981 


"1962 
"1942 
+1920 
+1898 


"1875 


*1850 
*1825 
*1799 
‘1771 


1743 


“1714 
+1684 
"1653 
*1622 


"1589 


+1203 
*1167 


‘1131 ; 


*1095 
*1060 


+1026 
"0992 
*09g00 
"0929 
*ogoo 


4293 
5202 
3402 
8922 
1803 


2096 
9865 
5185 
8144 
8841 
7387 
3908 
8540 
1432 
2748 
2664 
1372 
9076 
5997 
2368 


2 4284 


2006 
8844 


1974 
6331 
2191 
327 


6327 


Ay (k) 
00063 6588 
*O127 2985 
*OIgo0 
"0254 
"0317 


39099 
4440 


gIt5 


*O35I 2 
"0444 | 
"0507 
"0570 
70033 


-0606 
0758 4030 
‘0820 
"0882 
"0944 


*1005 
*1006 


3812 
3805 
*1127 0381 
*1187 3316 


“1247 2380 


*1306 7340 
"1365 7 
"1424 39 


/ 
*1482 5153 





*1540 1 


1924 


6972 


°1597 
"1053 
+1709 6081 
+1764 
“1819 


8954 
5280 
*1873 4735 
+1926 6981 
*1979 1660 
+2030 8395 
*2051 6797 

31 6436 
*2180 6864 


3 7596 


5 S10g 





| Agglk) 


| o-0000 3183 
| 0001 2734 
| +0002 8656 
| +0005 
| +0007 


0957 
9644 
| *OOII 4730 
| ‘O15 6228 
| +9020 4156 
5534 


933 


0025 
0031 


*0035 6730 
+0046 0603 
| "0054 1033 
"0062 8056 
"0072 1710 
"0082 2036 
godso 
2892 
3524 


1036 


"0092 
|} *Olo4 
*or16 
*O129 
| *0142 5485 
| *o156 
| *or7I 
| +0187 

"0203 


6950 
5493 
1197 
4147 


*0220 


238 


4433 
2154 
7415 
0276 0340 
1045 


0256 


0296 
| "0316 9665 
| +0338 6359 

0301 
| *0384 
; *0405 


4603 


6520 


*0433 
"0459 7040 
*0456 6143 
4845 
3494 


7201 


“0514 


| 0543 


2401 


2150 


0573 
+0604 
0636 3180 
6045 
"0704 1522 


0669 


"0740 O515 
*O777 4224 
| *0816 4252 
+0857 3225 
"0go0 6327 














Tables of the Bennett Functions A,,,,,(k) 



























































































150 Agg(k) Aj3(k) | Agg(R) 
0°02 "0423 77 0°0021 2111 00000 3183 0°0000 0011 | 00000 0000 
04 "0421 8705 "0042 3649 - *0000 0085 *0000 C000 
06 "0418 702 "0063 4043 *0000 0287 | *O000 0002 
08 0414 0084 2723 +0000 0680 *0000 0007 
10 "0408 64 O104 QI121 *0000 1329 | *0000 0017 
12 401 i 5 2675 ‘OOIr 4178 “0000 2298 | -0000 0035 
14 039 . O145 "OO15 5205 "0000 3653 | "0000 0064 
16 384 64 “0164 "0020 2409 *0000 5459 *0000 OIIO 
18 0374 } “0184 "0025 5732 *0000 7782 | *0000 0177 
20 0363 1 "0202 "0031 5107 - ‘0001 C691 | +0000 0270 
22 ‘Cc 544 0220 8462 "0035 0459 "OOO! 4253 | 0000 0397 
24 0337 6344 "0235 3464 "0045 1705 — *OOOI 5537 "0000 0565 
26 5524 "0052 8754 — *0002 3614 | *0000 0781 
28 55 ‘O 3190 *OO61 1505 *0002 9556 *0000 1057 
30 "0293 0453 "0286 6935 — -0069 9849 | — "0003 6435 *0000 1400 
| 
32 0276 7679 "0301 2640 "0079 3666 | "0004 4328 | ‘0000 1823 
34 "0259 9135 "0314 9552 "0089 2829 "0005 3310 *0000 2338 
36 0242 5060 "0327 5135 *0099 7200 *0006 3461 "E000 2959 
38 0224 8131 "0339 7080 ‘O110 6629 "0007 4862 *0000 3700 
40 0206 7458 *0350 6282 *O122 0960 *0008 7595 *0000 4578 
2 0158 4583 *0360 5369 *O134 0020 ‘0010 1748 | *0000 5610 
44 > 0471 0369 3988 - +0146 3630 *OOII 7409 | ‘0000 6817 
46 OSI 6116 "0377 1813 *O159 1595 "0013 4072 *0000 5219 
48 "0133 2527 "0353 8543 - °O172 3711 *OO15 3632 | ‘0000 9840 
50 “OII5 0732 "0359 3904 *O155 9757 *OOI7 4359 “OOO! 1705 
52 9097 17° "0393 76054 *O199 9501 *OO1g 7O50 “OOO! 3844 
54 0079 065 "0396 958: "0214 2093 "0022 1723 ‘OOo! 6286 
56 9062 6514 0398 9520 *0228 go7I 0024 8524 | ‘ooo! 9067 
58 0046 2308 "0399 7325 "0243 8354 "0027 7574 "0002 2224 
60 030 50 399 2901 "0259 0241 "0030 9003 | ‘0002 5798 
62 oo1s 5881 "0274 4413 "0034 2948 | *0002 9836 
64 OOO! 5 "0290 0531 "0037 9553 *0003 4390 
66 I 462 0305 8231 | 0041 8976 | +0003 9518 
68 23 4004 *O321 7122 "0046 1383 | ‘0004 5284 
79 7 "0337 6788 "0050 6956 "0005 1762 
72 3 70 0353 6777 0055 5892 | *0005 9034 
74 I go52 "0369 6604 ‘0060 8406 | +0006 7196 
76 0058 724 *0385 5743 "0066 4736 | *0007 6355 
78 064 11¢ *O401 3621 "0072 5145 | ‘0008 6639 
80 068 0257 *0416 9610 0078 9928 | 0009 8193 
82 14 0316 2215 "0432 3017 “0085 9417 ‘OOII 1192 
84 0302 4270 "0447 3072 "0093 3997 | ‘0012 5844 
86 070 8198 “0461 8908 ‘OIOI 4112 *QO14 2403 
88 2068 837 0475 9539 "0110 0289 | +0016 1185 
90 65 5052 0257 0688 0489 3817 *O119 3170 “0018 2588 
| 

92 60 9498 "O241 1135 "0502 0377 *O129 3555 *0020 7142 
94 0055 35 0225 0735 0513 7533 "0140 2495 | 0023 5575 
96 0048 997 “0209 2369 0524 3088 | "0152 1455 | 0026 8962 
"98 "0042 2915 0194 OIIO "0533 3865 "0105 2735 "0030 QII2 
I-00 "0030 0253 *O150 1205 *0540 3790 | *O150 1265 "0036 0253 








CORRIGENDA 


SOUTHWELL, ‘SOME EXTENSIONS OF 
“RAYLEIGH’S PRINCIPLE”’ 


(Vol. VI, Pt. 3, 1953) 


. 259, section 3, second line following equation (7): ‘A,’ should be ‘A,’ 

. 261, section 6, equations (11): A? and A® in first integrals (on left) should be 
deleted 

. 264, section 9, equation (19)A: ‘2(1-+-e)’ should be ‘—2(1+.e)’ j 
section 10, third line following equation (25): for ‘nearer A?’ substitute ‘nearer 
the smaller of A? and A2.,’ 

. 265, section 11, fourth line following equations (26): ‘A,’ should be ‘A,’ 

. 266, section 13, first equation (29): ‘> (C?/A?)’ should be ‘ > (C3)/> (C2/Az)’ 


. 267, second line: ‘A;,’ should be ‘w, 


CORRIGENDA 


R. V. SOUTHWELL AND G. VAISEY, ‘ON SOME 

EIGENVALUE PROBLEMS OF EXCEPTIONAL 

DIFFICULTY, EXEMPLIFIED BY A CASE OF 
ELASTIC INSTABILITY’ 


(Vol. VI, Pt. 4, 1953) 


. 453, second line from end of page: ‘in general’ should be ‘sometimes’ 

. 460, last line should read: ‘negative quantity, thus decreasing |,,| and increasing 
A_,|’ 

. 462, line 3, right-hand side of equation (23): ‘wy,’ should be ‘A.wy’ 

equation (25): ‘(cf—bg)*y”’ should be ‘(cf— bg)y?’ 

. 464, section 15, equations (24)A: ‘A?’ should be deleted in last 3 lines 

. 465, section 17, second line from end of page: ‘overestimates’ should be ‘may over- 
estimate’ 

. 467, section 19, second line: ‘will make’ should be ‘will be needed to make’ 
Fig. 5: legend on left should be ‘pp (force at P for wp = 1)’ 

. 469, line 7: ‘fq’ should be ‘@p’ 

. 471, section 26, second line: ‘section 12’ should be ‘section 10’ 

. 473, first line: ‘unit forces’ should be ‘unit force’ 

. 474, section 30, fifth linet ‘0-1’ should be ‘0-01’ 

line following equation (43): ‘(37)’ should be ‘(37)A’ 


w? 


. 475, third line: ‘section 7’ should be ‘section 6’ 


. 479, twelfth line: ‘section 88’ should be ‘sections 88 and 92a’ 
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